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Abstract 


A new method for orbit prediction, which is as accurate as numerical methods and as fast 
as analytical methods, in terms of computational time, is desirable. This paper presents 
Kolmogorov Arnol’d Moser (KAM) torus orbit prediction using Simplified General Per¬ 
turbations 4 (SGP4) and Two-Line Element (TLE) data. First, a periodic orbit and its 
Floquet solution are calculated. After that, perturbations, which are on the order of 10 5 
and smaller, are added to the periodic orbit plus Floquet solution. Then, the low eccentric¬ 
ity KAM torus is least squares fitted to the SGP4 and TEE data. The performance of the 
theory is presented in various ways. The new method is approximately five times more ac¬ 
curate for the best fits and three times more accurate for mean fits comparing to SGP4 and 
TLE. History of TLEs and KAM torus theory can be used to make accurate orbit predic¬ 
tions, which is conceptually similar to extrapolation. In addition, the new method may rival 
numerical methods and it can be used for collision avoidance calculations, and formation 
flight applications. However, high eccentricity, polar and critical inclination, air drag, and 
resonance problems should be addressed. 
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KAM TORUS ORBIT PREDICTION FROM TWO LINE ELEMENT SETS 


I. Introduction 

Artificial satellites are indispensable parts of the daily lives of people. Both civilian and 
military world depend on satellites orbiting Earth at very high speeds. There are many ap¬ 
plications of satellites such as communication, navigation, Earth observation, and weather. 
The position of a satellite should be known within a certain limit to establish communi¬ 
cation between ground stations and the satellite. The accurate positions of satellites and 
debris are important for collision avoidance purposes. Orbit determination methods yield 
the position of a satellite at a given time in the future. 

1.1 Motivation and Background 

Bright moving objects in the sky have intrigued mankind since ancient history began. 
Ancient people worshiped heavenly bodies as Gods. They developed timekeeping systems 
to facilitate their daily lives. They used stars to navigate in deserts and seas. They made 
observations to figure out the motion of celestial bodies. It was Johannes Kepler who 
first developed the three laws of planetary motion in 1619. Kepler’s empirical findings, 
based on Tycho Brahe’s planetary observations, were explained by Isaac Newton in his 
book “Principia” in 1687. Isaac Newton used his three laws of motion to derive Kepler’s 
laws of planetary motion. Kepler’s equations serve as a complete solution to the two- 
body problem (2BP), which determines the mutual gravitational interaction of two bodies 
under the assumption that there is no third body affecting them. The solution to the 2BP 
assumes that there is no perturbing force affecting the two body motion. The formulation 
of the three-body problem (3BP) was provided by Joseph Louis Lagrange in 1772. The 
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3BP is a special case of the n-body problem (NBP). It is the problem of calculating mutual 
gravitational interaction of three bodies. 

At the end of the 18 th century, remarkable advances occured in perturbation theory due 
to the logarithms invented by John Napier in 1614 and the efforts in modeling the Earth’s 
gravitational field. Perturbations can be defined as the divergence from the normal motion. 
The motion of a body under the influence of small accelerations diverges from the two 
body motion. There are three different methods to find perturbations to the two body mo¬ 
tion: special perturbation techniques, general perturbation techniques, and semianalytical 
techniques. Special perturbation techniques rely on numerical integration techiques and re¬ 
quire high computational power. Numerical integration methods provide orbit predictions 
with a level of accuracy in meters. General perturbation techniques make use of analyt¬ 
ical solutions to the perturbation problem and do not require much computational power. 
However, predictions made by general perturbation techniques are not as accurate as that of 
numerical methods. Semianalytical techniques combine special and general perturbation 
techniques to compensate for their disadvantages [60]. 

Precise orbit prediction is of great importance in the space age. Accurate positions 
of operational satellites help us to communicate with them. Antennas are useless if they 
are not pointing to the right direction, unless they are omnidirectional. Nonoperational 
satellites and debris pose threats to operational satellites. Accurate orbit estimation is very 
important for collision avoidance because wrong maneuver decisions cannot only result 
in the loss of millions of dollars but also create thousands of pieces of debris. Moreover, 
accurate predictions for the reentry of returning space objects are very important for two 
reasons. First, space objects returning from space might pose threats to people. Second, 
space assets returning from space might contain hardware or software important to the 
national security. Therefore, U.S. Space Surveillance Network (SSN) detects, tracks, and 
catalogs debris, operational and nonoperational satellites every day. SSN has more than 


2 



23,000 objects orbiting Earth in its database. [54]. 

North American Aerospace Defense Command (NORAD) publishes Two-Line Ele¬ 
ment (TLE) sets of some satellites to promote situational awareness in space. Special 
General Perturbations 4 (SGP4) is the current method of NORAD to propagate orbits from 
raw observational data coming from radars. However, SGP4 is accurate in the vicinity of 
epoch time. SGP4 is based on general perturbation techniques. Therefore, it suffers from 
errors that emerge from omission of higher order terms in the equations of motions, and its 
predictions are accurate only in kilometers. SGP4 requires observational data to make cor¬ 
rections once every few days. However, it is very fast in terms of computation time. Special 
perturbation techniques are used if more accurate predictions are needed. For example, the 
success of a collision avoidance intrinsically depends on accuracy in the positions of ob¬ 
jects in danger of collision. The downsides of numerical integration methods are truncation 
errors in the long run due to the limited word length of computers, high computational cost, 
and difficulty in determining whether the resultant orbit is the one that is desired. More¬ 
over, the breakup of satellite 1961-Omicron in 1961 proved that heavily dependence on 
numerical methods that requires high computational power for the sake of accuracy is not 
reliable. In 1961, the breakup tripled the number of debris to be tracked. The Naval Ordi¬ 
nance Research Calculator (NORC), which was at the Naval Weapon’s Laboratory (NWL) 
in Dahlgren, was insufficient to cope with the surveillance data flow and the orbital up¬ 
dates. Recently acquired IBM 7090, which was three times faster than NORC in terms of 
computational power, and incredible human effort processed the data in a few days. That 
catastrophic event has been second to none since then. However, lessons were well learned. 
It is crucial to have accurate orbit determination models that don’t require as much compu¬ 
tational power as special perturbation techniques [26]. Therefore, a new method for orbit 
prediction, which is as accurate as numerical methods and as fast as analytical methods, in 
terms of computational time, would be very desirable. 
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1.2 Problem Statement 


This work will answer the question of whether best low eccentricity Kolmogorov Arnold 
Moser (KAM) tori fitted to the orbits, using TLE histories of satellites, produce more accu¬ 
rate predictions for low eccentricity orbits than SGP4 does. Moreover, this effort presents 
an approximate accuracy analysis for the low KAM torus orbit prediction because both 
SGP4 and TLE suffer from intrinsic inaccuracies. The actual accuracy that can be achieved 
by this theory can be analyzed using raw observational data. 

A previous study by Frey showed that it is possible to extract KAM torus basis fre¬ 
quencies from SGP4 and TLE sets. Frey concluded that low eccentricities of Hubble Space 
Telescope (HST) and Thor Rocket Body caused difficulties in the construction process of 
KAM torus, and the difficulties in the TLE curve fitting for the two Delta Rocket Bod¬ 
ies are due to air drag. Frey used a modified Laskar frequency algorithm, developed by 
Wiesel, to determine KAM torus basis frequencies [28]. Laskar fequency analysis is used 
to determine the stability of orbits in dynamical systems where the energy is conserved. 

Because the modified Laskar frequency algorithm cannot model low eccentricity orbits 
to the desired degree of accuracy, Wiesel developed a new theory for low eccentricity Earth 
satellite motion to construct KAM torus. This work will make use of this new theory to 
build KAM torus for low eccentricity orbits. This KAM torus construction method includes 
geopotential, order and degree 20, and air drag perturbations [69]. 

1.3 Approach 

This effort uses three different main programs to prove KAM torus orbit prediction 
from TLE sets provides more accurate predictions in the long run than SGP4 and TLE sets. 
The first main program is a public domain SGP4 program written in C++ by Vallado [61]. 
It outputs predicted positions and velocities for an orbit period and some orbital elements 
by making use of TLE sets of satellites. TLE sets of satellites with low eccentricities are 
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obtained from www.space-track.org once every 3 to 4 days for a period of 2 months. NO¬ 
RAD imports corrections to SGP4 predictions with observational data coming from radars 
once every 3 to 4 days. The second main program, developed by Wiesel, builds a low ec¬ 
centricity KAM torus including zonal, sectoral, and tesseral geopotential perturbations, and 
air drag perturbations. The third main program is a least squares fitting program that fits 
a low eccentricity KAM torus to SGP4 and TLE data by iteratively correcting modal vari¬ 
ables, which are KAM torus coordinates and their linearizations, and basis frequencies of 
the KAM torus, and subsequently outputs residual mean square values. The methodology 
is based on the concept, which is very similar to that of the least squares method, that it is 
possible to obtain more accurate results with more data, and SGP4 is more accurate when 
nearer to the epoch time. The author has developed some scripts to analyze all low eccen¬ 
tricity satellites, which have eccentricities smaller than 10 2 , and inclinations not close to 
critical and polar inclination. There are 7,938 satellites and pieces of debris correspond 
with the above description, which can be publicly obtained from www.space-track.org. 
However, the author has pseudo-randomly chosen 1500 of these 7,398 near-Earth objects 
as his test cases. This effort serves as a prelude to converting the TLE sets catalog from 
SGP4 to low eccentricity KAM torus theory. 

1.4 Results 

This effort proves that the low eccentricity KAM theory is a better substitute for SGP4. 
The new theory is five times more accurate for the best fits, and three times better for mean 
fits. Most of the low eccentricity orbits of the near-Earth objects can be modeled by the low 
eccentricity KAM torus method. However, resonance, higher eccentricity, polar and critical 
inclinations, and air drag issues should be addressed. The results are certainly promising. 
Specifically, the new method can be used for collision avoidance calculations, and forma¬ 
tion flight applications. The new method provides with a set of numerical algorithms that 
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may rival the numerical methods in accuracy and the analytic methods in computational 
speed. This work also presents a rigorous performance analysis of the theory, which will 
be hopefully used as a reference for future studies. 

1.5 Overview 

This document is organized in five chapters. Chapter II will present previous studies 
conducted by Wiesel and his masters and PhD students on KAM torus orbit determination. 
It will also give an overview of the history of analytical orbit modeling in the United States 
space surveillance system, satellite dynamics formulated in Hamiltonian mechanics, KAM 
theory, the geopotential and air drag as perturbing forces, and least squares fitting. Chapter 
III will discuss the methodology and the data to be analyzed. It will also outline predicting 
positions and velocities from TLE and SGP4, building low eccentricity KAM tori, and 
fitting the low eccentricity KAM tori to the orbits obtained from SGP4 and TLE. Chapter IV 
will present the results and limitations of KAM torus orbit prediction method. It will also 
compare the accuracy of the new method to that of SGP4 and TLE. A general summary of 
the work, conclusions, and recommendations for future studies will be presented in Chapter 
V. 
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II. Background 


This chapter provides background information on current orbit determination methods, 
near-Earth satellite motion, and KAM theory. First, a brief history of analytical orbit mod¬ 
eling in the United States Space Surveillance System is presented. Then, the importance of 
the space situational awareness is presented. Next, Hamiltonian mechanics, transformation 
theory, orbit perturbations, and KAM theory are described. Then, the application of KAM 
theory to celestial mechanics is presented. After that, the previous studies of Wiesel and 
his masters and PhD students on the applications of KAM theory to Earth satellites are 
provided in chronological order. Finally, a summary of this chapter will be given. 

2.1 A Brief History of Analytical Orbit Modeling 

After the launch of Sputnik in 1957, United States started to track space objects. Today, 
SSN database has more than 23,000 space objects orbiting earth [54]. Earth orbiting objects 
are cataloged by analytical orbit models. The analytical methods have been mostly devel¬ 
oped by the Air Force Space Command and the Naval Space Command. The development 
and improvement of orbital models and algorithms span from 1957 to today [26]. 

The necessity of knowing orbital positions of space objects emerged due to military 
concerns in 1957. The Air Force used it not to confuse a missile with an object orbiting 
Earth, and the Navy used it to warn the fleets against space reconnaissance. The potential 
benefits from artificial satellites led to great interest in the accurate positions of the satellites 
for the civilian world. Therefore, the catalog for Earth orbiting objects was created, and it 
is still in use today [26]. 

The first formal catalog was created at the National Space Surveillance Control Center 
(NSSCC). The control center was located in Bedford, Massachusetts. The observational 
data were provided by 150 different sites. There were 4 different types of observations, 


7 



which were categorized by their source and content. Table 1 shows details of different 
observational data types used from the year 1957 to 1963 [26]. 

Table 1. Observational Data Types [26] 


Observation Type 

Content 

Source 

Type 1 

2 angles and slant range 

Radars 

Type 2 

2 angles 

Baker-Nunn cameras, 
telescopes, binoculars, 
visual sightings 

Type 3 

Azimuth 

Direction finders 

Type 4 

Time of closest approach (Doppler) 

Radars, radio receivers 
(for transmitting satellites) 


Observations were processed by an IBM-709 computer to update orbital elements in 
the catalog, which were some type of mean orbital elements. Then, updated orbital data 
were processed by another program to yield three different products, which were used to 
make observations for the next time by ground sites. The first product, called bulletin, was 
the precursor of the TLE sets. The most important components of a bulletin were longi¬ 
tude, time of ascending pass, and revolution number of a satellite. These were generated 
approximately a week in advance. The second product was used to calculate all satellite 
passes, which could be seen from an observing site. The third product was very similar 
to the second product. However, the third one was used by the U.S. Navy and U.S. Army 
observation fences, which were composed of a number of radars. Therefore, the third prod¬ 
uct was used to yield the intersection point of orbital plane of a satellite and vertical plane 
spanned by radar beam instead of some look angles and ranges [26]. 

In 1961, U.S. Navy constructed the Naval Space Surveillance System (NAVSPASUR), 
which detected and cataloged Earth orbiting objects mostly without human intervention. 
The processing unit for NAVSPASUR was located at Dahlgren, Virginia. The system was 
a continuous wave multi-static interferometer, which observed LEO satellites 4 to 6 times a 
day. There were 3 transmitters and 6 receivers, which spanned from San Diego, California 




to Savannah, Georgia as 9 individual sites located in the Southern States. The system 
was called “Fence”. The command of the Fence, which was responsible for 40% of all 
observations made by the Air Force, was passed to U.S. Air Force in 2004. The U.S. Air 
Force shut down the system in October 2013. A new fence is planned to be operational in 
2017 in the Marshall Islands [19]. 

Special perturbation techniques were thought to be used in the Fence until it was dis¬ 
covered that it took the NORC at the NWL several hours to update one orbit. In 1961, 
the catalog was transferred to an IBM-7090 computer, which reduced the time required 
to update one orbit to approximately 1 minute. That was very promising for the imple¬ 
mentation of highly numerical methods, which yield more accurate predictions. However, 
the breakup of 1961-Omicron tripled the number of objects to be tracked. The number 
of detectable objects exceeded the capability of the NORC. The IBM-7090 and incredible 
human effort processed the data in a few days. Although special perturbation techniques 
provide with more accurate predictions, the beakup of satellite 1961-Omicron proved that 
accurate analytical models are needed in case computers fail to process the data [26]. 

In 1959, Dirk Brouwer and Yoshide Kozai provided two different solutions to the same 
problem, which was Earth satellite motion under the influence of the zonal harmonics J2, 
J3, J4, and J5 [9, 38]. Modern orbit prediction methods used in the U.S. Space Surveillance 
System include fundementals of the solutions provided by Brouwer and Kozai. In 1961, 
Brouwer and Gen-Ichiro Hori added atmospheric drag effect to the 1959 Brouwer solution 
[10]. However, the atmospheric model was computationally heavy for the computers of 
the time due to the series expansions in the scale height. In 1960, a group of people under 
the project SPACETRACK developed an atmospheric density modeling, which solved the 
series expansion problem for artificial satellites. Max Lane developed an analytic orbit 
model based on the atmospheric density modeling of the SPACETRACK in 1965 [40]. 
Lane and Cranford improved the analytic orbit model by implementing analytic density 
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model instead of empirical one in 1969 [42]. In 1963, Lyddane contributed to the analytic 
orbit modeling by solving small divisors of the eccentricity and sine of the inclination 
by defining the perturbation theory in terms of Poincare variables instead of Delaunay 
variables [49]. 

Conversion from theory to practice for analytical orbit models took place in Dahlgren, 
Virginia, and Colorado Springs, Colorado. NAVSPASUR implemented the 1959 Brouwer 
solution with Lyddane’s contribution, which was known as PPT (Position and Partials as a 
function of Time), on an IBM-7090 computer in 1964. PPT adopted semi-empirical drag 
model developed by Richard H. Smith because the atmospheric model of Brouwer and Hori 
was computationally heavy for the computers of the time, and it was too early for the Lane’s 
analytic orbit model. The change in the eccentricity could be solved using Equation 1 and 
Equation 2. 


e D = e a (l ~e 0 2 ) — 

do 


( 1 ) 


4 a 0 d C) 

a ° = ~ 3 ^ 2 ) 


(2) 


Numerically, mean motion term in PPT is similar to Kozai’s mean motion, which has the 
zonal secular perturbation rate of mean anomaly drived by Brouwer. 

NSSCC was moved to Colorado Springs, and named as Space Detection and Tracking 
System (SPADATS). SPADATS is the other location for the implementation of the ana¬ 
lytic orbital theory, which was called SGP (Simplified General Perturbations), developed 
by Brouwer and Kozai. Arsenault, Chaffee, and Kuhlman avoided small divisors of the 
eccentricity and sine of the inclination in SGP by defining the solution in non-singular 
terms and keeping only the most important ones in the model. They included long and 
short-period terms without the eccentricity effect from the Brouwer solution and the con- 
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vention that relates mean motion to semi-major axis from the Kozai solution [33]. The air 
drag effect used in SGP was very similar to Smith’s semi-empirical model. In 1969, Lane 
and Cranford developed an analytic atmospheric model for artificial satellites [42]. How¬ 
ever, computers of the time were inadequate to cope with calculations required for the new 
model. Therefore, the most important terms of the atmospheric model were implemented 
to the analytical orbit model, which is now called SGP4. Today, some flavor of SGP4 is 
used by NORAD [26]. 

Solar and Lunar gravitation, and resonance in the tesseral harmonics of Earth became 
important when the first satellite with a period of 12 hours was launched in 1965. Bruce 
Bowman developed a model which includes third body effects, and resonance in the tesseral 
harmonics of Earth in 1967 [7]. In 1977, Dick Hujsak developed another model which in¬ 
cludes everything from Bowman’s work. Moreover, his first-order model could be applied 
to geosynchronous satellites [32]. Hujsak’s model was combined with the SGP4 model, 
which is still in use today. In 1997, third body effects, and resonance in the tesseral har¬ 
monics of Earth were adopted from SGP4 by the Naval Space Command. The new method, 
which is still in use today, is called PPT3 [26]. 

2.2 Space Situational Awareness 

In less than six years from the launch of Sputnik 1 to 1963, 616 man-made objects 
accumulated in low Earth orbit (LEO). Those man-made objects were 76 payloads, 35 
rocket bodies, and 227 mission related debris. The U.S. Ablestar upper stage explosion 
created 184 pieces of objects which were half of all the objects cataloged by SSN in 1963, 
and 60% of all debris from this explosion is still orbiting Earth [54]. Today, 23,000 objects 
are cataloged by SSN and the majority of these objects are not operational satellites but 
debris and nonoperational satellites, which is approximately 95% of total number of objects 
being cataloged [59]. 
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Figure 1. Cataloged objects by SSN in 1963 (left) and in 2013 (right) [54] 

The rate of growth in the number of debri objects in less than six years was noticed 
by Ernest Peterkin, Head of the Systems Section of the Operational Research Branch at 
the U.S. Naval Research Laboratory. He was one of the first to calculate the number of 
objects per unit volume in LEO in Lebruary 1963. He derived population growth functions 
for the objects in LEO using estimates for launch rates, lifetimes of satellites, and satellite 
collisions. Peterkin’s conclusions showed that there would be 16,500 objects in LEO at 
the beginning of 2013. His predictions are very close to the number of objects cataloged 
by SSN on 1 January 2013. Peterkin concluded that space surveillance systems would 
have difficulty coping with the rapid growth in the number of objects orbiting Earth in 
the future. He stated that it is important to improve detection and data processing methods 
used in space surveillance systems to manage the increase in the number of objects orbiting 
Earth [54]. This effort is trying to find an accurate and a fast orbit prediction method 
which can help space surveillance systems to cope with the proliferation of objects orbiting 
Earth. Low eccentricity KAM theory yields more accurate orbit predictions than SGP4. In 
addition, this work shows that it is possible to extract valuable data, which are used by the 
low eccentricity KAM torus orbit prediction method to make corrections to low eccentricity 
KAM tori, from inaccurate TLEs in the long run. 

The Joint Space Operations Center (JSpOC), located at Vandenberg, is responsible for 
detecting, tracking, and cataloging all man-made objects orbiting Earth. JSpOC uses SSN 
which has 30 space surveillance radars and optical telescopes all over the world. The radar 
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Figure 2. Worldwide network of space surveillance sensors of SSN [55] 

sensors are dedicated to make observations for near-Earth objects, which are below 6000 
kilometers altitude. The electro-optic sensors are dedicated to make observations for deep- 
space objects, which are above 6000 kilometers altitude. Azimuth, elevation, range, and 
range rate are typical observation types obtained from a radar site. Azimuth and elevation 
or right ascension and declination are observations provided by an optical sensor. Observa¬ 
tions made by SSN include time tags [37]. Figure 2 shows locations of worldwide network 
of 30 space surveillance sensors. Observational data obtained by these sensors are fitted to 
the trajectories of orbiting objects, and the satellite catalog is continuously updated. SSN 
makes 380,000 to 420,000 observations every day. SSN implements a predictive surveil¬ 
lance method because all objects orbiting Earth cannot be covered by existing sensors, 
which requires efficient analytic orbit models. These observations are important for opera¬ 
tional satellites because collision analyses are made using orbit predictions in the satellite 
catalog [59]. 

This effort is beneficial for satellite communities which use NORAD TLE for propa¬ 
gating orbits of the objects that pose threat to their space assets because low eccentricity 
KAM theory can yield position and velocity data to create TLEs before NORAD publishes 
them. NORAD publishes a new element set when the difference of positions predicted by 
the current and the new element sets exceed 5 km with a 90% confidence interval. LEO 
satellites require more frequent updates due to the atmospheric drag. Maneuverable satel- 
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lites also need more frequent updates because orbital maneuvers cannot be modeled by 
SGP4 [37]. This work uses Vallado’s public domain SGP4 with differential corrections to 
determine the accuracy of the orbit predictions yielded by low eccentricity KAM theory 
because JSpOC version of SGP4 has no technical details and codes available to public. 

2.3 Orbit Prediction using SGP4 and TLE 

TLEs are some type of mean elements which are averaged, and smoothly changing over 
time or angle. Some period of time, mean anomaly, true anomaly, and eccentric anomaly 
can be used to average the elements [60]. TLEs are used by NORAD, and the majority of 
the satellite community. NORAD TLE is the default input for most of the commercialized 
satellite ground antenna control system and orbit analysis programs [46]. TLE is gener¬ 
ated from observational data flow from SSN using one of five orbit prediction methods 
of NORAD. The first one is SGP which was developed by Hilton and Kuhlman [31], and 
simplified by the work of Kozai in terms of gravitational and air drag effects [38]. SGP 
is used for satellites with a period less than 225 minutes. The second one is SGP4 which 
was developed by Cranford in 1970 [41], and simplified by adopting analytic atmospheric 
model of Lane and Cranford instead of empirical atmospheric model in 1969 [42]. Table 2 
shows parameters required to initialize SGP4. All of the parameters but mean motion are 
mean orbital elements defined by Brouwer [9]. Brouwer’s mean elements include only the 
effects of the zonal harmonics which are J2, J3, J4, and J5. SGP4 includes orbital mean 
elements and their linearizations defined in series expansions. Many assumptions are made 
for SGP4. For example, Earth’s gravitational potential includes only a few zonal terms, 
the atmospheric model is a static model with an assumption of exponential decay, and the 
third-body mass and resonance effects are partially added to the model [62]. Mean motion 
in the Table 2 includes only short-period oscillations similar to Kozai’s mean motion [38]. 
Long-periodic oscillations are masked by air drag, which gradually increases mean motion 
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Figure 3. TLE Example for NOAA-6 Weather Satellite 


over time, for LEO satellites [28]. 

Table 2. Parameters to Initialize SGP4 Propagation 


Symbols 

Description 

to 

Epoch time 

n a 

Mean motion at epoch 


Eccentricity at epoch 

io 

Inclination at epoch 

(0 0 

Argument of perigee at epoch 

o 0 

Right ascension of the ascending node at epoch 

M 0 

Mean anomaly at epoch 

B* 

Atmospheric drag coefficient 


The third one is SDP4 which is developed for orbits with a period more than 225 min¬ 
utes by Hujsak in 1979 [32]. SDP4 includes third-body, sectoral, and tesseral effects, 
whereas SGP4 doesn’t. The fourth one is SGP8 which is applicable to near-Earth satellites. 
SGP8 depends on the model of Lane and Cranford for gravitational and air drag effects 
[42]. The integration of differential equations are treated in a different way for SGP8 than 
SGP4. The fifth one is SDP8 which is a deep-space model. SGP8/SDP8 were introduced to 
decrease the effects of the deficiencies of SGP4/SDP4 associated with reentry and orbital 
decay [62]. A TLE should be used with one of these models because periodic variations 
are removed in a way that these five models can compensate for them. Figure 3 shows an 
example of a NORAD TLE with descriptions of each elements. 

Orbit prediction methods are divided into three different catagories: special perturba- 
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tion techniques, general perturbation techniques and semianalytical techniques. Special 
perturbation techniques use numerical methods to make predictions. However, numerical 
methods are prone to errors due to limited word length of computers. Moreover, long time 
steps for integrations result in round-off errors. Special perturbation techniques cannot be 
generalized to other orbits because the predictions are special to the orbit that is integrated. 
Special perturbation techniques provide predictions accurate in meters. However, they re¬ 
quire high computational power. All perturbations are calculated for every point in time 
in special perturbation techniques. General perturbation techniques use closed form solu¬ 
tions to equations of motions of orbits. The complexity of equations of motions requires 
omission of higher order terms in the equations. The accuracy of the predictions degrades 
rapidly due to these omissions. In addition, orbit predictions made by general perturbation 
techniques are not as accurate as those provided by special perturbation techniques. How¬ 
ever, general perturbation techniques provide orbit predictions without calculating every 
intermediate points to find the desired point in the future and take an average of all pertur¬ 
bations over time as one parameter which can be integrated into the equations. Therefore, 
they are very fast in terms of computation time. SGP4 orbit determination method is such a 
general perturbation technique. Semianalytical techniques combine advantages of analyti¬ 
cal methods and numerical methods. 

Equation 3 shows the relationship between the SGP4 model and the state vector of a 
satellite. 

y(t) = fsGP4(x 0 ,B*,t), (3) 

where /sgpa is the SGP4 model, y(t) is the velocity and position vectors at time t, x 0 is 
the mean orbital elements, which are shown in the Table 2 excluding t a and B*, at the 
epoch time. The epoch time is represented as t 0 [46]. Equation 4 shows how the ballistic 
coefficient, B , and the atmospheric reference reference density p G are related to NORAD 
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SGP4 air drag term B * [52]. 


B* = l -B Po (4) 

Equation 5 can also be used to produce the NORAD TLE given osculating elements and 
a good approximation for the NORAD TLE. Equation 5 shows the relationship between os¬ 
culating elements and the NORAD mean elements. B* is missing in the Equation 5 because 
SGP4 can’t model the air drag term. Therefore, B* needs to be calculated separately [46]. 
B* is important for the accuracy of the NORAD TLE. The omission of B value of 10 4 
could worsen the accuracy of SGP4 twice as much within approximately 3 days [45]. 

yi= fsGP4i(xi,...,x 6 ), (5) 

where i=( 1 ,• • • ,6), and y,- are osculating elements, x- L are the NORAD mean elements. Equa¬ 
tion 5 can be expanded around xf-, which are good approximations for the NORAD mean 
elements, using Taylor series expansion. Equation 6 and Equation 7 show the linear system 
of equations after taylor series expansion [45]. 

Ay/ = MAxj, (6) 

dfsGP4 l d/SGP4 l df S gP4 l 

dx^ < 7*2 

M= \ (7) 

dfsGP4 6 dfsCP4 6 dfsGP4 6 

dx^ <7c“ dxfc 

where Ay,- = y / — yf, and At, = x, — x“. Partial derivatives in matrix M are hard to solve due 
to the coupled nature of the SGP4 model. Newton’s quotient can be used to find the partial 
derivatives [25]. Equation 8 shows one of the partial derivatives in the difference quotient 
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notation [45]. 


dfsGPA, = /sGP4« + Axf ,X%, • • • ,jcg) - /SGP4(*1,*2,'' ' X) 
r7.v" Ax“ ’ 1 J 

where Ar" can be some small increment of percentage of jt". The NORAD mean elements 
can be approximated using a process called differential corrections. The iterative process 
is stopped when the difference between the given osculating elements and derived ones, 
which are (y,- — yf), are small enough for convergence. Equation 8 shows the iterative 
process to approximate the NORAD mean elements [45]. 

^ +1 =*?+EVfc-tf)’ (9) 

7=1 

where i=(l,-- - ,6). 

There are inevitable deficiencies in the SGP4 theory. Because SGP4 is an analytic or¬ 
bit model, it ignores higher order terms in the equations of motion. Moreover, sectoral 
and tesseral gravity field perturbations are not included in the model. The SGP4 model 
includes only the effects of the zonal harmonics which are J2, J3, J4,and J5. Therefore, the 
SGP4 model yields position errors of 2 km at epoch time [52]. The velocity predictions 
made by the NORAD TLE and the SGP4 are less accurate than the position predictions be¬ 
cause the rates of change in the orbital parameters, which are not accurate due to inevitable 
deficiencies in the SGP4 theory, are used to calculate velocities [28]. 

The accuracy of the Vallado’s SGP4 with differential corrections for near-Earth objects 
are on the order of magnitude of 100 km. The prediction accuracy worsens with decreasing 
altitude, and increasing eccentricity and solar flux. Table 3 and Table 4 show the SGP4 
orbit prediction accuracy of the near-Earth objects in terms of the root mean square values 
with respect to two different solar flux values [24]. 
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Table 3. The Orbit Prediction Accuracy of the Near-Earth Objects (E 10 = 100)[24] 


Altitude(km) 

Test Cases 

lday(km) 

3days(km) 

7days(km) 

15days(km) 

h<400 

57 

4 

10 

60 

300 

400<h<600 

168 

3 

10 

50 

100 

600<h<1200 

267 

3 

10 

20 

50 

1200<h<7000 

90 

2 

10 

10 

20 


Table 4. The Orbit Prediction Accuracy of the Near-Earth Objects (E 10 = 200)[24] 


Altitude(km) 

Test Cases 

lday(km) 

3days(km) 

7days(km) 

15days(km) 

h<400 

57 

10 

40 

300 

1000 

400<h<600 

168 

7 

30 

200 

400 

600<h<1200 

267 

6 

15 

70 

100 

1200<h<7000 

90 

2 

10 

10 

20 


2.4 Reference Frames 

Reference frames are important for classical and analytic mechanics. All equations 
of motions are written in some frame of reference. The first step in solving a dynamical 
problem is defining the reference frame. There are two different types of reference frames, 
which are inertial and rotating frames. There is no solution for any dynamical problem 
in classical mechanics unless an inertial frame is defined. This effort formulate satellite 
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Table 5. Reference Frames 


Acronym 

Reference 

Frame 

Type 

Origin 

I st axis 

2 nd axis 

3 rd axis 

RAN 

Radial, 

Along- 

Track,and 

Normal 

Rotating 

frame 

Space 

object 

CoM 

Toward 

radial 

direction 

(U) 

Normal 

to orbital 

frame 

(W) 

Along 

velocity 

direction 

(V) 

ECI 

Earth- 

Centered- 

Inertial 

frame 

Inertial 

frame 

Earth 

CoM 

Toward 

vernal 

equinox 

(7) 

Completes 
right- 
handed 
frame (£) 

Earth’s 

axis of 
rotation 
(0 

ECEF 

Earth- 

Centered- 

Earth- 

Fixed 

frame 

Rotating 

frame 

Earth 

CoM 

Toward 
prime 
meridian 
in equa¬ 
torial 
plane (x) 

Completes 
right- 
handed 
frame (y) 

Earth’s 

axis of 
rotation 

(z) 

ECNF 

Earth- 

Centered- 

Node- 

Fixed 

frame 

Rotating 

frame 

Earth 

CoM 

Toward 
ascend¬ 
ing node 
in equa¬ 
torial 
plane (£) 

Completes 

right- 

handed 

frame 

(V) 

Earth’s 

axis of 
rotation 
(0 


motion using 4 different reference frames, which are defined in Table 5 and Figure 4. 

Reference frames in Table 5 are helpful in different ways for this effort. RAN reference 
frame is used to make residuals between SGP4 and low eccentricity KAM theory predic¬ 
tions meaningful because it is hard to analyze the residuals in ECI frame. ECNF reference 
frame is important because orbits are not periodic in ECI frame although they are in ECNF 
frame. Periodic orbits are not periodic in ECI frame because orbital plane regresses due to 
the zonal potential. Moreover, the Keplerian frequency is defined in a frame which is tied 
to the node of the orbit. Thus, it can simply be generalized to ECI and ECEF frame because 
the zonal potential is symmetric about the earth’s polar axis. Greenwich referenced ECEF 
frame simplifies calculations because the gravitational field potential is constant in ECEF 
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frame. 


2.5 Analytical Dynamics 

Equations of motions are formulated using momentum and force in Newtonian me¬ 
chanics. Objects and forces applied on them are considered separately. For Newtonian 
mechanics, constraint forces, which are forces that do no work, may be needed to solve 
dynamical problems. However, analytical dynamics considers a dynamical system as a 
whole not separately. For analytical dynamics, equation of motions are formulated using 
kinetic energy and work, which are both scalar quantities. Therefore, constraint forces are 
not needed to solve dynamical problems. Dynamical systems in analytical dynamics are 
defined in generalized coordinates, which are finite, continuous, and differentiable with 
respect to time. In addition, generalized coordinates aren’t restricted to be physical quan¬ 
tities. Number of generalized coordinates for a system equals to the number of degrees of 
freedom that the system has [51]. Many complex dynamical systems can be solved without 
much effort using analytical dynamics. Figure 5 shows the evolution of transition from 
Newtonian to Fagrangian mechanics. 

Principle of Virtual Work 

I 

D' Alembert Principle 

i 

Hamilton's Law 

I 

Hamilton's Principle (Hamilton's Extended Principle) 


Lagrange's Equations 

Figure 5. Evolution of the Transition from Newtonian to Lagrangian Mechanics 

Principle of virtual work is based on static equilibrium of a system. It is the first varia¬ 
tional principle in mechanics [51]. Constraint forces are perpendicular to virtual displace- 
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ments. However, formulation of equations of motions relies on forces. Equation 10 and 
Equation 11 show the virtual work of the entire system equals to zero. 

_ N N 

5 W = £ • 8n + £ Fj ■ 8 n = 0, (10) 

i= 1 i= I 

N 

£ ■ 8rt = 0, (11) 

1=1 

where Fj are external forces, and F[ are constraint forces applied on the system. Sr, are 
virtual displacements. 

D’ Alembert’s principle is an extension of principle of virtual work to dynamics. D’ 
Alembert’s principle treats dynamical problems as if they were statical. Equation 12 shows 
the most general formulation of dynamics [51]. 

N 

£ (Fj — Pi) 8rj = 0, (12) 

i=l 

where pi is rate of change of the momentum, which is also called inertia force. 

Hamilton’s law is important because it depends on kinetic and potential energy instead 
of force and momentum. Both principal of virtual work and D’ Alembert’s principle are 
similar to Newtonian mechanics in formulating equations of motion because they all de¬ 
pend on force and momentum. The first step in transition from D’ Alembert’s principle 
to Hamilton’s law is integrating D’ Alembert’s principal equation over a finite period of 
time, and expressing momentum in terms of kinetic energy. The derivation of Hamilton’s 
law introduces Lagrangian, which is the difference between kinetic and potential energy. 
Equation 13 represents Hamilton’s law of varying action. Hamilton’s law of varying action 
is valid for non-linear, non-conservative, non-stationary dynamic systems. Variations of 
the generalized coordinates Sqj aren’t necessarily zero because the system state qj isn’t 
always known at time t\ and h, which shows the system is non-stationary. Variation of any 
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constant value equals to zero [39]. 



( 13 ) 


L = T — V, (14) 

where L is Lagrangian, T is kinetic energy,and V is potential energy. 

Hamilton’s principle is an extension of Hamilton’s law. For Hamilton’s principle, the 
generalized coordinates dqj are assumed to be known at end points t\ and * 2 - Hamilton’s 
principle is applicable when there are no non-conservative forces, whereas extended Hamil¬ 
ton’s principle is applicable when there are non-conservative forces in the system. Equa¬ 
tion 15 shows Hamilton’s principle, and Equation 16 shows extended Hamilton’s principle 
[39]. Hamilton’s principle is a variational principle which reduces a dynamical problem 
to a scalar integral independent of coordinates, which describes Lagrangian. Equations of 
motion are obtained by figuring out the conditions which make the scalar integral station¬ 
ary. Figure 6 represents true path and varied path, which the true path coincide with two 
end points t\ and h- The varied path is formed when the virtual displacements is applied, 
which by definition, has no change in 5r, [51]. 



0 


(15) 


8Ldt + / " 8Wdt = 0 


(16) 


Lagrange equations are n number of coupled second-order differential equations, which 
has n number of generalized coordinates. Lagrange equations are derived from Hamilton’s 
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Figure 6. True Path and Varied Path 


law of varying action as [39]: 


rt 2 _ N 

8L{q i ,q il t)dt+ 8W-Y{^8qj] 
Jn “l ch h 


= 0 


(17) 


After calculating the variation of the Lagrangian, 


di, dL „ dL 


—— A.dT 


£/ ( t 7 : s 4 i +-^- i 8q i )dt+ 


'h dqi dq { 


i d 4j 


= 0 


(18) 


After integrating Lagrangian terms by parts, and cancelling the trailing term, which shows 
the system could be non-stationary, 


y f 2 \ c L(—)- — -Q i 

“j Jq [dl dqi dqi 


8q/dt — 0, 


(19) 


where Q t is generalized force. The integrand can be equated to zero because the system 
can be considered stationary. Equation 21 shows Lagrange’s equations when there is no 
non-conservative force in the system because virtual work equals to zero. 


d ^ dL ^ dL 
dt dqi 7 dqi 


0 


( 20 ) 
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0 


( 21 ) 


d ^ dL ^ dL 
dt dip dqi 


2.6 Hamilton’s Equations 


The Hamiltonian formulation is an alternative to the Lagrangian formulation [16]. They 
both have the same physics. The Lagrangian formulation has n number of second-order 
differential equations, whereas the Hamiltonian formulation has 2n number of first-order 
differential equations. The Lagrange’s equations describe mechanics in terms of general¬ 
ized coordinates and velocities, whereas the Hamilton’s equations rely on the generalized 
coordinates, qp and momenta, pt [29]. Equation 22 shows the equation for the generalized 
momenta. The generalized momenta can be described as linear functions of the general¬ 
ized velocities. Conversely, the generalized velocities can be shown to be linear functions 
of the generalized momenta [51]. The Legendre transform of the Lagrangian yields the 
Hamiltonian. Equation 23 shows the Hamiltonian. 


dL(qi,qj,t ) 

dqt 


( 22 ) 


H (Piqit) = Y<<iiPi-L(q,qd), (23) 

i 

where H(p,q,t) is Hamiltonian, which can fully describe the motion [51]. 

The Hamilton’s formalism are more powerful than Lagrange formalism. The apparent 
advantage of the Hamilton’s equations over Lagrange equations is that time derivatives of 
the variables are on the left side of the equations, which may help to figure out first integrals 
of motion. Moreover, the Hamilton’s equations are favorable in representing the solution 
in geometrically because 2n-dimensional representation is advantageous in representing 
not only one path but totality of paths, which are all possible solutions. In Hamiltonian 
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phase space, the velocity of every point can be determined uniquely, whereas the motion 
can start in any direction, and paths may intersect in the Lagrangian configuration space, 
which makes totality of paths confusing [51]. The canonical equations of Hamiltonian can 
be derived from the Hamiltonian as [29]: 


dtf _ dH , dH J 

dH — > — dqj+ > -=r— dp t + -z—dt 

i uC[i • opt dt 


(24) 


After substituting Equation 23 into Equation 24, 


dH = l^qidpi - Y, ^ dc li ~ ~^dt 


(25) 


and from Equation 22 it follows as, 


dL 


dH = Yd idpi~YPi dc li ~ ~^dt 


(26) 


After equating each term in Equation 26 to each term in Equation 24, 


dH 


dH 


li j Pi — -j 
opi dcfa 


(27) 


Equations 27 are the canonical equations of Hamilton, which are 2n number of first-order 
differential equations. 


2.7 Canonical Transformations 

Hamiltonian formalism doesn’t introduce an efficient calculation tool by itself. The 
freedom in choosing the physical quantities as coordinates and momenta is beneficial for 
solving difficult dynamical problems. Changing coordinates helps to reveal much about any 
dynamical system. In Hamiltonian mechanics, the Hamiltonian is constant if it doesn’t ex- 
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plicitly depend on time, and the momenta are constant if the associated coordinates, which 
are called ignorable or cyclic coordinates, are missing in the Hamiltonian. However, it is 
necessary to start from the beginning when the change of coordinates is required. More¬ 
over, for Hamiltonian formalism, the coordinates can be changed to anything as long as 
they describe the system, but the momenta are dictated by the Lagrangian, see Equation 22. 
One method for changing coordinates without starting from the beginning by conserving 
the structure of Hamiltonian dynamics is canonical transformation [64]. 

Assume Equation 28 represents new canonical coordinates, and Equation 29 represents 
new momenta: 

Qi = Qi(qi,Pi,t) (28) 


Pi = Piiqi.pi.t) 


(29) 


Also assume a new Hamiltonian K = K(Qi,Pi,t), and Hamilton’s equations as: 


dK . _ dK 

dPi 1 Pi ~~dQ 


The old variables q t and p t , and the new variables Qj and P t are supposed to describe the 
same dynamical system. The Hamilton’s principle can be used to show that the variation 
of both integrals equal to zero as: 

, N 

8 / (^Pi4i- H )d t = 0, (31) 

J i= 1 


r N 

8 ('£PiQi~ H )dt = o, (32) 

2 7=1 

where L can be deduced from Legendre transform as L = (Lfli Pi4i ~H). Although the 
variations of both integrals equal to zero, the integrands aren’t equal. The total time deriva- 
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tive of an arbitrary function F separates two integrands from each other as: 


8 I £ Pi4i - H (. Pi , qut) ~ Y, P iQi 
dF 

+K(Qi,Pi,t) - — ]dt = 0 


( 33 ) 


The variation vanishes at end times as: 


5 f 2 ^dt = 8{F{t 2 )-F{t x )) = 0, (34) 

Jt\ dt 

where F is the generating function, which is a function of independent 2n variables. Ta¬ 
ble 6 shows the four possible forms of the generating function F. The generating function 
F needs to be known in order to apply the transformation. F can be obtained from the 
relationship between old and new coordinates, and Table 6, see Wiesel [64]. 

Table 6. Canonical Transform Relations 


Fi(q,Q,t) 

p 2 (q,P,t) 

p 3 (p,Q,t) 

Hp-Pd) 

Pl ~ d qi 

n ■ — FF 

Pl ~ dq t 

dF$ 

di - ~Jfd 

_ _ dp4 

h ~~ dpi 

p. — dF] 

F ~ dQj 

O — dFi 
'd' ~ rJP, 

p — dFi 
p ~ dQj 

n. — dF* 
dl — dpj 


The new Hamiltonian K and the old Hamiltonian H are related by Equation 35 as: 

dF] 

K(Q,P,t)=H(q,p,t) + ^~ (35) 

2.8 Hamilton-Jacobi Theory 

The easier solution to a dynamical system can be accomplished when the more coor¬ 
dinates are missing in the Hamiltonian. Moreover, any momenta missing from the new 
Hamiltonian dictates that its conjugate coordinate is constant. If the new momenta, K, 
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equals to zero, all the new coordinates and momenta are constant [64]. Equation 36 shows 
that the new Hamiltonian equals to zero, which is called the Hamilton-Jacobi equation, and 
Equation 37 shows that the new coordinates and momenta are constant [29]. 


H{quPi,t) + 


dF 

~dt 


0 


(36) 


Qi 


Pi 


dK 

dPi 

dK 


0 

0 


(37) 


It is favorable to use a generating function which is a function of the old coordinates, 
qi, and the new momenta, P,. The generating function, F 2 , from Table 6 is such a function. 
Therefore, Equation 36 becomes [29]: 


H 


dF 2 

QU 3 j t 


dt 


(38) 


Equation 38 is known to be Hamilton’s principal function. The solution to this equation is 
denoted by S as [29]: 

S = S(qi • • • < 7/15 CL\ • • • , ( 39 ) 

where a, are n number of constants of the integration because S is the solution to the first- 
order differential equation. Therefore, one of the variables of the solution has to be an 
arbitrary constant added to the solution S. Using the analogy of the physical description 
of the generating function, the constants of the integration, oq, can be chosen as the new 
momenta, Pj [29]: 

Pi = oc, (40) 
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After applying transformation for the fo generating function from Table 6 [29]: 


dS{qi,ai,t) 

P ‘ dqi 

Continuing with the second half of the transformation [29]: 

„ _ « _ dS(qi,CCi,t ) 

Mi — Pi — ^ 


( 41 ) 


(42) 


Both 0Ci and /3, can be obtained by evaluating partial derivatives with the known initial 
conditions. The Hamilton-Jacobi equation can be solved as [29]: 


q = q((Xi,pi,t) 


(43) 


The Hamilton’s principal function helps to transform the old variables to new constant 
coordinates and momenta. A solution to the dynamical system is obtained when solving the 
Hamilton-Jacobi equation. The relationship between the Hamilton’s principal function and 
the Hamilton’s principle can be shown by examining total time derivative of the Hamilton’s 
principle function, S [29]: 


dS y-, dS dS 
dt jdqi^'+dt 

After inserting Equation 36 and Equation 41, Equation 44 becomes: 


(44) 


=Y,Piqi~ H = L (45) 

Therefore, a constant differs the Hamiltonian principle function from Hamilton’s principle 
as: 

S = I Ldt + constant (46) 

For Hamilton’s principle, the solution to a dynamical system is obtained by solving the 
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definite integral of Lagrangian, L. Equation 46 represents the same integral in indefinite 
form which can be used to solve a dynamical system. 

Hamilton’s principle function can be separated in two parts when the Hamiltonian 
doesn’t explicitly depend on time as [29] 


S{qu (X{,t ) = W (<?/,«/) - oqf, (47) 

where VE(^/,cq) is the first part of the solution, which depends on the old coordinates 
and the new momenta, and appears only when the Hamiltonian is constant. W{qi, ccf) is 
called Hamilton’s characteristic function. Equation 47 is the solution to the Hamilton- 
Jacobi equation [29]: 

f + «te.f) = 0 (48) 

After substituting Equation 47, the Equation 48 becomes [29]: 

dW 

H(c li —) = a u (49) 

dc/i 

where one of the constants of integration, cq, equals to the constant value of the Hamilto¬ 
nian, H. The new Hamiltonian, K, equals to oq because H doesn’t depend on time [29]: 

K = a i (50) 

After calculating Hamilton’s equations for the new Hamiltonian, K [29]: 

P, = CCj 

i= 1 (51) 

i ^ 1 


Pi = 
Qi = 


dK n 


dH 

r)a, 


= 1 


= 0 
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After denoting the generating function by W (, qpPi ), the solution becomes [29]: 


Qi 


Qi 


t + Pi = 


dw 

d OC\ 



i ^ 1 


( 52 ) 


where Q\ is not a constant of the motion, and time , t, is a coordinate and the Hamiltonian 
is its conjugate momenta [29]. 


2.9 Action-Angle Variables 


The periodic motion is of great importance in physics. The frequencies of the motion is 
mostly more desirable than the other properties of the orbit. A variation of the Hamilton- 
Jacobi process is used to obtain the frequency of the periodic motion. The action variables, 
J, are chosen instead of the new momenta of the Hamilton-Jacobi equations. /, are inde¬ 
pendent functions of the a, of the Section 2.8. Figure 7 shows an example of the two types 
of periodic motion. For oscillation, the system repeats its track for every point because q 
and p return to its original values after one period. For rotation, the position coordinate q 
is an unbounded angle of rotation, which grows indefinitely with a period of q 0 , whereas q 
is bounded for oscillation [29]. 

Ji can be defined as [29]: 

Ji = j pidqi (53) 

After applying the canonical transformation, Equation 53 becomes [29]: 


Ji 


dW(qj, eg) 
dqi 


dqii 


(54) 


where the solution to the Equation 54 is J\ as independent functions of a,. Substituting J, 
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Figure 7. The Orbit of Phase Space for Periodic Motion [29] 


for a, yields the characteristic function W as [29]: 


W = W(qi---q n ;Ji---J n 


(55) 


Also, the angle variables, ft)/, are the conjugate coordinates to the J, [29]: 

dW 


Pi = 


ft), = 


dqi 

dW 

~dJi 


(56) 


The Hamiltonian is a function of only i, because w, are cyclic [29]: 


H = H(J l --J n ) = a l 


(57) 


The Hamilton equations of motions for the new variables become [29]: 


dH n 
Ji = — 3 — = ° 

<9ft)/ 

dH 


(58) 
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where v, are constant functions of The solution becomes [29]: 


/, = constant 


(Oi = Vit + A 


( 59 ) 


where the solution doesn’t introduce any advantage over the a, coordinates from the Sec¬ 
tion 2.8. The advantage of this process is that the Action-Angle variables transformation 
yields the frequencies, V/, of the periodic motion without finding a complete solution. If a 
system is periodic, the energy of the system can be defined in terms of /, by Equation 53. 
Next, the frequencies, V/, of the system can be obtained by Equation 58. The to, are conve¬ 
niently called angle variables due to Equation 59 [29]. 

The fact that the frequency v, are related to qt can be proved by investigating the change 
in one of the angle variables, CO,-, when one of the coordinates, q,, completes its one cycle 
of oscillation or rotation. The Equation 60 shows that the change in the angle variable is 
caused by an incremental increase in one of the coordinates, qj [29]: 


A (Ot = 



(60) 


where 8 C0j is the infinitesimal change due to the incremental increase in one of the coordi¬ 
nates, qj [29]: 

8(Oi—^r—^-dqj (61) 

dqj 

Also, by Equation 56 and Equation 61, Aft), can be written as [29]: 


A ft), 


d ft), 
dqj 


dqj — 


d 2 W , 

dq^J lj 


(62) 


The /, in the denominator can be taken out of the integral because 7, are independent func- 
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tions of only a, [29]: 


. s rsw . s 

,/.7 / 


By Equation 53, Aft), becomes [29]: 


(63) 


. Sj, 

AMi = Hi = l 

= 8ij i ^ j 


(64) 


If r,■ is the period of the motion related to qi. Equation 65 becomes [29]: 


Aft), - 1 = V/Tj 
1 

V/ = - 

T; 


(65) 


2.10 General Canonical Transformations 

It is not an easy task to create a generating function for every canonical transformations. 
A different approach to canonical transformations can be initiated as [64]: 


x T = {qi, Pi] 


( 66 ) 


Next, the Hamilton’s equations of motion becomes [64]: 


Pi = 


c)H 

dpi 

dH 

dqi 


(67) 
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Equation 67 can be defined in terms of phase space state vector by Equation 66 [64]: 


x — Z 


Z = 


dH 

dx 

0 / 
-1 0 


( 68 ) 


where Z satisfies the properties [64]: 


-z = z l =z 


(69) 


Nothing new has been introduced except representing the Hamilton’s equations in terms 
of the phase space state vector. The transformation of the phase space from the old vari¬ 
ables, x, to the new variables, y, can be accomplished as [64]: 

X = f(y) (70) 

Then, the new Hamiltonian, K , becomes [64]: 


K(y) = H(fW), 


(71) 


where the old Hamiltonian, H , doesn’t explicitly depend on time. The transformation, /, is 
valid as long as it is a canonical transformation. Therefore, the conditions which make / a 
canonical transformation should be examined. The new coordinates and momenta should 
comply with the Hamilton’s equations as [64]: 


y — Z 


dK 

dy 


(72) 


Converting Equation 72 to Equation 68 may reveal these conditions. Equation 73 shows 
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the time derivative of the transformation, and Equation 74 represents the gradient of the new 
Hamiltonian [64]: 


i = 


df. 

dy y 


(73) 


c)K d . (df y dH 


dy dy H ^ f ^ \ dy) dx 


(74) 


Combining these equations with Equation 72 yields [64]: 


. ,3/ l. df r dH 

y = { d~ y ) x=z( Ty ] at 


(75) 


Then, Equation 75 yields x as [64]: 


• t d f>7t S f>TdH 

x = ( Ty )z( Ty ] at 


(76) 


Comparing Equation 68 to Equation 76 immediately yields [64]: 


(—L)z(—) t = z 

[ dy )Z[ dy ) Z ’ 


(77) 


where the first-order partial derivative is called the Jacobian matrix J. The J matrix is 
symplectic when it satisfies [64]: 

JZJ T = Z (78) 

Therefore, the transformation, /, is a canonical transformation if its Jacobian matrix is 
symplectic. 


2.11 Orbit Perturbations for Near-Earth Satellites 

Orbit perturbations are the deviations from the two-body orbit motion. Two-body mo¬ 
tion assumes that secondary body orbits around a point of mass or spherically symmetric 
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sphere [15]. However, Earth is neither a point of mass nor spherically symmetric. The 
equatorial bulge, continental blocks, ocean basis, and mountain ranges result in gravita¬ 
tional field deviations from the two-body point of mass. Another important perturbation 
for near-Earth satellites is air drag. Air drag effect is bigger when a satellite has smaller 
mass with bigger cross section. Air drag acceleration is a non-conservative force, which 
determines the orbit lifetime of a satellite in Low Earth Orbit (LEO) [64]. Table 7 repre¬ 
sents three categories for Earth orbits. LEO satellites are primarily affected by air drag and 
non-sphericity of Earth, and both MEO and GEO satellites are mainly perturbed by solar- 
radiation pressure and third-body effects [60]. This work includes air drag, and sectoral and 
tesseral harmonics perturbation accelerations. The effects of zonal harmonics perturbation 
are included in the periodic orbit, see Section 3.2. 

Table 7. The Three Categories of Earth Orbits [60] 


Orbit 

Altitude (km) 

Low Earth Orbit (LEO) 

h < 800 

Medium Earth Orbit (MEO) 

800 < h < 30,000 

Geostationary Orbit (GEO) 

h = 35,780 


The gravitational field equation, which is also called Poisson’s equation, is the funda¬ 
mental partial differential equations for gravitational fields [64]: 


V 2 V =4nGp 1 (79) 

where G is the gravitational constant, p is mass distribution, V 2 is Laplacian operator, and 
V is gravitational potential function. Equation 79 defines the gravitational potential of a 
sphere. The gravitational potential can be calculated by Equation 79 when mass distri¬ 
bution, p, is provided. However, the gravitational potential outside the sphere is needed. 
Equation 80 represents the infinite series of the gravitational potential function outside a 
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Zonal Harmonics 
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Figure 8. The Visual Depictions of Spherical Harmonics [60] 

sphere, which has a radius of 2?© [64]. 

72=0/72=0 ® (80) 

x ( C nm cosm(j) + S nm sinmty) 

where /l is the gravitational parameter, n and m are the order and degree of the expansion, 
respectively, 7?© is the radius of Earth, P™ are the associated Legendre functions, the 0 
is the geocentric latitude, (f> is the longitude, and C nm and S nm are gravitational constants, 
which can be obtained from the gravitation models. Both C nm and S nm define the shape of 
the gravitational field [64]. All spherical harmonics in the gravitational potential field of 
Earth, which are zonal, tesseral, and sectoral harmonics, can be represented by Equation 80 
[60]. Figure 8 shows the visual depictions of three types of spherical harmonics. 

Newtonian point of mass potential can be obtained from Equation 80 when m,n = 0. 
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The first term is free of longitude and latitude terms because cos m(j) = 0, sin inO = 0, 
Coo = 1, and the associate Legendre function Pq = 1. 

V = --, (81) 

r 

where V is the Newtonian point of mass potential. 

The geopotential expansion yields zonal harmonics when n is the order and m — 0 
is the degree. The smallest zonal harmonic, which has n = 1 order and m = 0 degree, 
is considered to be zero because it shifts the center of mass of Earth parallel to North- 
South axis. Any frame which doesn’t have its origin at the center of mass is not practical. 
Therefore, the first zonal harmonic equals to zero, J\ = 0, when the origin is defined at the 
center of Earth. The second zonal harmonic specifies the oblateness of Earth at the equator. 
For the Earth’s potential expansion, the J 2 is second to the Newtonian point mass potential 
in strength. The Newtonian potential is 10 3 times bigger than the the 72. The contribution 
of J 2 term to the gravitational potential can be represented as [64]: 

V 2 o = ^ 1 (3cos 2 e-\) 1 (82) 

where J 2 — —C 20 — 0.001082. For the second zonal harmonic, J 2 , mass concentration is 
bigger at the equator than the poles. The higher order zonal harmonics add more terms to 
the potential, which model the irregularities of the mass distribution along the longitudes. 
See Figure 8 for the visual depiction of zonal harmonics [64]. 

The sectoral harmonics are the terms with equal orders and degrees in the potential 
expansion, n = m. They are functions of only longitude, and they change signs across 
the longitudes. The first sectoral harmonic, n = m = 1, changes the location of center of 
mass away from the origin of a reference frame to another one in the equatorial plane. 
Therefore, Sn = Cn =0 because it is always desirable to have center of mass at the origin 
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of a reference frame. The second sectoral harmonic, n — m = 2, contributes to the potential 
as [64]: 

V 22 °= C 22 COS 2 (j) + S 22 sin 2 (j ), (83) 

where (j) is the longitude. The second sectoral harmonics, C 22 and S 22 , are on the order of 
10 6 . J 2 is 10 3 times bigger than C 22 and ^ 22 - The mass concentrations of Earth at Eurasia 
and the Americas are bigger than Pacific and Atlantic ocean basis due to the second sectoral 
harmonics. See Figure 8 for the visual depiction of sectoral harmonics [64]. 

The tesseral harmonics are the case when n / m / 0 in the geopotential expansion. They 
create a sphere with rectangular square-tiled boards on it. The number of bands in latitude 
is n — m+ 1 because of P™ (cos 0) dependence. The terms , C nm cos m(j) and S nm sin m(j), 
disappear for 2m meridians. Therefore, there are 2m bands in longitude. See Figure 8 for 
the visual depiction of tesseral harmonics [64, 60]. 

The gravity models, such as EGM96, include the coefficients, C nm and S nm . The gravity 
models are built from terrestrial or satellite measurements. The satellite measurements are 
orbit dependent. Therefore, the accuracy of a gravity model depends on the number of 
satellites with different orbits. This effort uses the EGM96 gravity model. The EGM96 is 
developed by the University of Texas at Austin, the Defence Mapping Agency, Ohio State 
University, and the Goddard Spaceflight Center. The EGM96 includes 30 different satellite 
measurements and terrestrial gravity measurements. It is a complete model with degree and 
order 360 [60]. However, this work uses only order and degree 20 because it provides high 
accuracy with low computational power requirement. Figure 9 represents the variations in 
the geopotential. 

The air drag is the other perturbation included by this work. The air molecules behave 
individually, which is called the free molecular flow regime, at the altitude of satellites. 
The air molecules impact the motion of the satellites. Equation 84 represents the drag law 
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Figure 9. GRACE Gravity Model [53] 


Temperature Variations in the Atmosphere 
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Figure 10. Temperature Variations in the Atmosphere 


[64]: 

a d = -\— P VY, (84) 

2 m 

where C d is the drag coefficient, which is 2.2 for flat plane models, and between 2.0 and 
2.1 for spheres for the most part, p is the atmospheric density, A is the cross-sectional area 
of the satellite, V is the velocity of the satellite relative to the air molecules, and m is the 
mass of the satellite. The drag coefficient, C d , specifies the vulnerability of a satellite to 
the air drag, and it is dimensionless. The atmospheric density at the altitude of a satellite 
is represented by p, which is difficult to calculate due to the interaction between Earth’s 
upper atmosphere and Sun. The cross-sectional area, A, is also difficult to calculate because 
the attitude of the satellite must be known. Moreover, it is almost impossible to determine 
the attitude of a tumbling satellite [60]. 

Solar flux and geomagnetic storms heat up the thermosphere. Therefore, the upper at- 
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mosphere expands, which increases the number of air molecules. During Solar maximum 
the emission of Extreme Ultraviolet radiation (EUV) increases because Solar flare and 
Coronal Mass Ejection (CME) events increase during Solar maximum. Charged particles 
coming from the night side of the Earth’s atmosphere during geomagnetic storm interact 
with air molecules, and increase their energy, which expands the upper atmosphere. Ge¬ 
omagnetic storms cause delayed expansion in the atmosphere because storms first hit the 
day side of the Earth’s atmosphere. Both Solar flux and geomagnetic storms increase the 
atmospheric density in the upper atmosphere. Figure 10, which is plotted using Mass Spec¬ 
trometer - Incoherent Scatter (MSIS-E-90) model data, shows the temperature change in 
the atmosphere due to the diurnal variations and the Solar cycle, and Figure 11, which is 
plotted using SGP4 and TLE prediction data, represents the impact of air drag on the Delta 
1 rocket body over 138 days, which has an altitude of 482 km. 


6865 

6860 
E 

6855 

t 6850 
E 

| 6845 

6840 

6835 

23000 23020 23040 23060 23080 23100 23120 23140 23160 
Julian date 


DELTA 1 R/B 



Figure 11. Air Drag Effect on the Delta 1 Rocket Body 

Although the atmospheric density is unpredictable, it is very important for the accuracy 
of the orbit predictions. Many atmospheric models, which are either static or time-varying, 
has been developed. Every model is either based on physical models, which are built 
by combining conservation laws and atmospheric-constituent models or developed from 
in-situ measurements and satellite-tracking data. Every atmospheric model is different in 
terms of speed, accuracy, and applicability. Therefore, there is no model which provides 
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best results for all applications. This effort uses the 1962 Standard Atmosphere. The 
routine for the atmosphere model requires the altitude of each layer, molecular weight, and 
molecular temperature. Moreover, the temperature within each layer changes linearly, see 
Regan and Anandarskarian [27]. It is an ideal, static atmospheric model at a latitude of 45° 
during moderate Solar activity [60]. 


2.12 Numerical Integration Methods 


There are many numerical integration methods which have been developed to solve 
ordinary differential equations. Most of these methods have been successfully applied in 
celestial mechanics. However, there is no numerical integration method which yields the 
best solution to every problem pertaining to the motion of satellites. The most important 
numerical integration methods for orbit computations are Runge-Kutta methods, multistep 
methods, and extrapolation methods. Runge-Kutta methods are easy and applicable to wide 
range of different problems. Multistep methods are very efficient, but they need to store 
previous data points. Extrapolation methods are highly accurate [56]. 

This effort uses 4 -order Hamming predictor-corrector method, which is basically a 
multistep numerical integrator. It requires four values of the state vector. However, only 
the initial conditions are know at the beginning. The other three values can be obtained by 
a process called Picard iteration. Equation 85 shows the Picard iterative process: 


y n+ i(x) =}’„+ / f(t,y n (t))dt 

Jx ° (85) 

y n 0 ) =y n 


If the Picard iteration converges, the four values of the state vector are extrapolated to 
give the next value for the state vector, which is the predictor part. Next, the extrapolated 
state is corrected using better values of the equations of motion, which is the corrector 
part. Although predictor-corrector methods require complex algorithms, they are favorable 
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because of their stability. 


2.13 KAM theory 

The French mathematician and physicist Henri Poincare found out that 3BP is unsolv- 
able because of small divisors problem. He suggested that N-Body Problem is mathemati¬ 
cally unsolvable. Birkhoff and Smale proved Poincare’s earlier findings, and they suggested 
that some nonlinear systems aren’t solvable. The primary question was: When does a sys¬ 
tem behave chaotically? The KAM theory, which was announced by A.N. Kolmogorov in 
1954, and proved by J. Moser and V.I. Amol’d in the early 1960’s, is the solution to this 
question from the Hamiltonian aspect. The KAM theory explains what happens when an 
integrable Hamiltonian is perturbed. It is developed to overcome the difficulties in pertur¬ 
bation theory about small divisors [50]. The KAM theory states that an integrable, which 
means n constants of the motion are known, Hamiltonian has a phase space motion which 
lies on a ^-dimensional torus in 2/7-dimensional phase space, where n is the number of in¬ 
dependent coordinates [47, 6]. The qua si-periodic motion can be represented by action and 
angle variables, see Section 2.9. The quasi-periodic orbits describe integrable motion on 
the invariant torus. If any trajectory is on an invariant torus in phase space, it will stay on 
the torus [5]. The quasi-periodic motions have n number of frequencies, which is dictated 
by the Hamilton-Jacobi theory [65]. 

The fundamental equation for the KAM theory can be shown as : 

H(l, 6) = H 0 (l) + 0), (86) 

where H is the perturbed Hamiltonian, H\ is the perturbing Hamiltonian, H 0 is the in¬ 
tegrable Hamiltonian, e is a small real number, e << 1, and I and 0 are action-angle 
variables. H 0 and H\ must be real analytic functions, which are infinitely differentiable 
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Figure 12. A Visual Depiction of 1-Torus [50] 


functions. The action variable I is constant and 0 is a linear function of time because H 0 
is only a function of I, see Equation 27. Figure 12 represents 1-torus, and Figure 13 shows 
2-torus where the action variables, I, are constant and the angle variables, 0, are functions 
of time. 

If the solutions for the perturbed Hamiltonian Equation 86 lie on //-dimensional tori, 
there are new action-angle variables as [57]: 


tf(I,0) =//'(!') 


(87) 


Equation 87 represents the generating function, S, for the Hamilton-Jacobi transformation, 
see Section 2.8 [57]: 


dS(l') 


de 


( 88 ) 



01 ' 


New Hamiltonian is a function of only I’ because the 0 variables are cyclic, and the 
Hamilton-Jacobi equation becomes [57]: 


H(—,e) = H\ I') 


(89) 
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Figure 13. A Visual Depiction of 2-Torus 

Power series expansion in £ can yield a solution to the generating function, S, as [57]: 

S = S 0 + eS l + £ 2 S 2 + ••• (90) 


Substituting Equation 89 in Equation 88 assuming S 0 — I'-O because Equation 87 yields 
1 = 1' and 9 = 9' after the assumption [57]: 


TT . . dS i 2 dS 2 

H °V + e lw +e Je 
+ £//l(I/ + e W 


+ ■■■) 

+ ---,9)=H’(l') 


(91) 


The expansion of Equation 91 for small £ keeping only the first-order terms yields [57]: 


H 0 (l') + £ 


dHo 

dt 


dS i 

~de 


+eh 1 ( i',e) 


H'( I') 


(92) 


Then, H i (If 9) and Si (I', 9) can be represented as Fourier series [57]: 


Hi = ■ 9) 

m 

Si = J^S l . m (l')exp(im ■ 9) 

m 


(93) 


where m is a vector with n rows consisting of integers. Substituting Equation 93 in Equa- 


47 






Invariant Tori 
(Regular Motion) 


Resonant Torus Distorts 
(Chaotic Motion at "X") 


Figure 14. Elliptic Islands and Hyperbolic Fixed Points 


tion 92 yields [57]: 



(94) 


where V 0 is the unperturbed ^-dimensional frequency vector by Equation 59. The infinite 
sum must converge to provide a solution. Moreover, the same condition must be satisfied 
for S 2 , 1 S 3 , ■ • ■ if £ is expanded to higher orders. Although a method of successive ap¬ 
proximations to S has been shown here, the proof of KAM theory depends on complex 
successive approximations which converge much faster. 

The convergence of the sum depends on the denominator, m ■ V 0 (I / ). This problem is 
called the small divisors problem. The action variables, I, specify the resonant tori, which 
behave chaotically when £ > 0 , for the unperturbed system when m ■ v 0 (I') = 0. On the 
other hand, the nonresonant tori must satisfy the Equation 94 [57]: 


|m- v| > A r (v)|m| 9 1+1 ) 


(95) 


where m are vectors consisting of integers, the absolute value of m, |m|, is the sum of all m 
vectors up to number n, and K( v) > 0 is an arbitrary number independent of m. Moreover, 
m can’t be the zero vector for nonresonant tori [57]. Therefore, it is highly likely that a 
perturbed, nearly integrable, and periodic system is described by an invariant tori in phase 
space [ 1 ]. 

Another way to describe the resonant tori, which helps to visualize it, is by winding 
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number, a. Figure 13 can help to visualize of the ratio of two frequencies, Vi and Vo- 
Assume Vi is the frequency seen from top view, and V 2 is the frequency seen from cross 
sectional view. If the winding number, a = V1/V2, is a rational number, the tori associated 
with this winding number are called resonant tori. The tori with irrational winding number, 
a , are called nonresonant tori. For nonresonant tori, “almost all” orbits are preserved 
because “the KAM theory restore a measure of continuity to chaos” [50]. Figure 14, which 
consists of two surface of section plots, represents what happens if a resonant torus breaks 
up when a perturbation is introduced to the system. According to the Birkhoff’s theorem, 
alternating elliptic and hyperbolic fixed points are formed when the resonant torus breaks 
down. The two stable elliptic islands, which is represented by O , and the two unstable 
hyperbolic points, which is represented by X, are shown in Figure 14. 

2.14 KAM Theory Applications 

After KAM theory was published, astronomers applied it to astronomical models be¬ 
cause the motions in the solar system are comparatively bounded. However, application 
of KAM theory to astronomical models didn’t yield promising results because of the lim¬ 
itation on the size of the perturbation parameter, e, which is the ratio of masses. In 1963, 
Arnold endeavored to establish the existence of KAM tori for NBP, which he was partially 
successful. The primary question was whether there were the initial position and velocities 
of the bodies that keep the distance of the bodies from each other bounded for all the time 
in the NBP [2]. A combination of KAM theory and computer-assisted techniques named 
interval arithmetic, which is a computational technique that controls the errors introduced 
by numerical computations in a special way, was successfully used to prove the existence 
of KAM tori for the Restricted, Planar, Circular 3-Body Problem (RPC3BP) by Celletti and 
others. In 1997, Celletti and Chierchia proved the existence of quasi-periodic tori with a 
frequency close to the average frequency of Ceres for the Sun-Jupiter-Ceres problem [12]. 
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Locatelli and Giorgilli proved the existence of KAM tori describing the secular motions of 
Jupiter and Saturn for the obserbed values of the parameters [48]. Celletti and Chierchia 
investigated a truncated RPC3BP model for Sun, Jupiter, and Asteroid 12 Victoria. The in¬ 
variant tori bounding the motion of Victoria was built successfully for the astronomic value 
of the Jupiter-Sun mass ratio [13]. For further information, see a brief history of KAM tori 
for NBP [14]. 

The numerical applications of KAM theory originated from Binney and Spergel’s pa¬ 
per which spectra obtained from Fourier series of the coordinates of numerically integrated 
orbits under the influence of galactic potentials can be expressed as sums of integer mul¬ 
tiple of fundamental frequencies. Binney and others realized that orbits of iV-dimensional 
galaxies defined in phase space can be represented as V-tori [4, 6]. McGill and Binney de¬ 
veloped a method to build tori for Hamiltonians including general gravitational potentials. 
This method is based on the idea of distorting the analytic tori of a toy Hamiltonian into the 
desired tori using generating functions [18, 6]. Binney and Kumar generalized the method 
for obtaining the frequencies and angular variables related to the tori least-squares fitted to 
any Hamiltonian [3, 6]. Kaasalainen and Binney further refined the torus-fitting process of 
the method in order to include the case of non-rotating bar, which admits two major orbit 
families instead of one [36, 6]. Kaasalainen and Binney introduced point transformations 
into the method in order to solve the problem of a toy potential being too different from 
its target potential [35, 6]. Then, Kaasalainen showed that the method can be applied to 
globally chaotic regions as well [34, 6]. 

The only work related to this effort in terms of the application of KAM theory for ob¬ 
jects orbiting Earth has been done by Wiesel and his masters and PhD students. Their 
efforts are represented in a chronological order. Wiesel showed that Earth satellite orbits 
under the full geopotential effect are likely to be KAM tori in a reference frame rotating 
with Earth. First, he determined the frequencies of the orbit using a Laskar frequency 
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algorithm. Next, he obtained Fourier series coefficients from numerically integrated tra¬ 
jectories. Then, He fitted numerically integrated orbit to the multiple Fourier series, which 
defines the torus, using least squares fitting technique. The spectral analysis of almost 
all orbits, excluding chaotic ones, yielded three basis frequencies. Moreover, Wiesel con¬ 
cluded that other dynamical effects can be added as perturbations because of the single 
point construction of the Hamilton-Jacobi separation of variables solution [65]. His next 
paper related to the application of KAM torus theory to Earth satellites is the one that 
compares two different KAM torus construction algorithms, which are the least squares 
fitting of a KAM torus to a numerically integrated orbit and the extraction of the Fourier 
coefficients of the KAM torus from numerically integrated Fourier transform. An efficient 
KAM torus construction method, which yields the torus that passes through the given initial 
conditions, is important for applications, such as formation flight of satellites, stationkeep¬ 
ing satellite constellations, and compressing ephemeris data for navigational satellites [66]. 
Derbis tried to apply the KAM theory to precise satellite data from the GPS satellites. She 
applied a Fast Fourier Transform (FFT) to obtain and identify discrete frequencies from 
orbit data. However, some orbits showed inconsistencies related to their basis frequencies. 
She concluded that the inconsistencies in the third frequency, ( 03 , which is the apsidal re¬ 
gression rate, caused by the resonance in orbits of the GPS satellites [22]. Little tried to 
build KAM tori from orbit data of the Gravity Recovery Climate Experiment (GRACE) 
and Jason-1 satellites. He applied a modified Laskar Fourier transform algorithm to obtain 
the basis frequencies of the orbits. Although the residuals of the KAM torus for Jason-1 
satellite were as low as 1 km over a 30 day period, the KAM torus construction process 
failed for GRACE satellite. Little concluded that air drag and errors related to Fourier 
coefficients resulted in the poor results for GRACE satellite. Craft investigated the appli¬ 
cability of the KAM theory for the satellite formation flight in the full geopotential with an 
order and degree of 20. The KAM torus method yielded promising results for the satellite 
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formation flight applications. Orbits of satellite clusters provided with drift rates in the 
nanometer to micrometer range over 60 day integration interval when they were separated 
by 10 to 100 meters [20]. Bordner tried to apply the KAM theory to the GPS satellite 
orbits. His research questions were whether the KAM theory can increase the accuracy 
of GPS satellites and whether it can reduce the burden on GPS operations. However, the 
methods for constructing KAM tori failed to yield desired accuracy for operational GPS or¬ 
bits. He suggested that improved methods are needed to deal with complexities related to 
the maneuvering operational satellites [6]. Wiesel’s third paper related to the KAM theory 
compared the KAM tori built from 2BP, SGP4 model, and numerically integrated orbit with 
a degree and order of 20. He showed that a KAM torus is similar to a full analytic perturba¬ 
tion theory in many ways. However, for a KAM torus, frequencies aren’t approximations 
of the perturbation series expansions, and the torus is built numerically. Moreover, the 
comparison of numerical integration and the KAM torus, which yielded residuals smaller 
than 4 m over 10 years, showed that trajectories in the full geopotential are KAM tori [68]. 
Yates investigated the motion near a reference torus in order to compensate for errors in the 
reference KAM torus due to dissipative forces, such as air drag, lunar/solar mass effects. 
He suggested that routine phase angle updates and stochastic offsets to the reference torus 
can improve the accuracy of the KAM torus perturbed by dissipative forces. He showed 
that most stochastic effects can be modeled to predict the motion near a reference torus. 
Moreover, the low eccentricity of International Space Station (ISS) is proved to be a big 
hurdle in estimating the torus coordinate offsets [70]. Hagen studied the effects of air drag 
and lunar mass perturbations on a reference KAM torus. He showed that the KAM theory 
is also applicable when air drag and lunar mass perturbations are included in the system. 
However, it is proved that an accurate Jacobian, which performs the linear transformation, 
is of great importance, and Jacobian is only valid when the reference torus coordinates 
are tightly aligned with the perturbed torus coordinates [30]. This effort includes air drag 
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as perturbations to the reference KAM torus. Frey studied feasibility of constructing the 
KAM torus using TLE and SGP4. Frey’s work, in a sense, proved that the KAM torus 
can be constructed from observational data. Two of his test cases, which are the two Delta 
rocket bodies, failed due to air drag and inaccuracies of the TLE data. The other two test 
cases, which are the Hubble Space Telescope (HST) and the Thor rocket body, showed 
that it is possible to extract accurate KAM torus basis frequencies although they suffered 
from their small eccentricities. Frey concluded that small eccentricities can’t be modeled 
with a modified Laskar frequency algorithm developed by Wiesel [28]. This effort uses a 
different method which builds the reference KAM torus for orbits with low eccentricities 
using periodic orbits and Floquet theory. 

This effort is based on Wiesel’s recent paper, “A Theory of Low Eccentricity Earth 
Satellite Motion”, and it builds on the results of previous efforts related to the application 
of KAM theory to Earth orbiting satellites. Wiesel considered the Earth satellite motion 
in terms of periodic orbits and Floquet theory. Periodic orbits in the zonal potential are 
known to be nearly circular, excluding the ones with critical inclination [58]. The main 
advantage of applying this method is that perturbations are in the order of 10 5 instead of 
10 3 because periodic orbits includes the effects of the Earth’s oblateness. The perturba¬ 
tions caused by sectoral and tesseral potential terms, air drag, and third body mass effects 
can be added to the periodic orbit. The current method is a hybrid of numerical and an¬ 
alytic methods. It has numerical sets of algorithms that may match numerical integration 
in their accuracy. Moreover, Once the theory files are constructed, which takes approxi¬ 
mately 1 minute per orbit, it has the the usual advantages of general perturbations because 
it provides the results at the time of interest without having to perform a long propagation 
[69], 
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2.15 Summary 


Space surveillance systems are under great stress due to the number of objects orbiting 
the Earth. Current launch rates and operational lifetimes of current satellites show that 
space surveillance systems need to cope with more objects soon. The increase in the spatial 
density of debris leads to an increase in the probability of collision. Figure 15 represents 
the growth of the cataloged satellite population during the past 15 years. Therefore, a new 
method which is as accurate as numerical methods and as fast as analytical methods is 
needed. The KAM theory for Earth orbiting objects yield promising results specifically 
for space debris because they are non-maneuverable,and lightly perturbed. Once a KAM 
torus is built for a non-maneuverable object, its motion will be bounded by the torus until 
a dissipative force affects its motion. Formation flight of satellites is another potential 
application of the KAM tori. Craft showed that tight formations yield secular drift rates 
between satellites from 4 nanometers to 1 micrometer per second [20]. Moreover, The 
KAM tori can reduce the operational burdens of GPS because the ephemerides provided by 
the orbital tori will have a much longer useful life [6]. Indeed, the KAM tori has numerous 
potential applications because they provide accurate predictions with long lasting validity 
for the time of interest as fast as an analytic method can provide. However, computationally 
expensive long integrations for building a torus, and perturbations that shift a reference 
torus to the nearby torus have been the main problems for the theory. Hopefully, this 
effort mostly overcomes these issues and the KAM torus orbit determination will serve as 
a prelude to converting the TLE catalog. 
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Figure 15. The Growth of the Cataloged Satellite Population over 15 Years [55] 
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III. Methodology 


This chapter describes the data, the method, and the analytic process used to perform 
this study. This effort is intended for small eccentric orbits, which were one of the main 
issues for the KAM torus construction process for the previous efforts. Although this effort 
has been built on the previous works, it introduces a completely new method for construct¬ 
ing the KAM torus, which was developed and coded by Wiesel. In addition, the new theory 
is a complete orbit determination method which retrieves the TLE data for two months from 
www. space-track. com servers, propagates each TLE for one orbit, build the KAM torus 
theory file using some of the SGP4 mean elements, and fits the low eccentricity KAM torus 
to the SGP4 and TLE predictions using least squares without human intervention. Although 
converting whole TLE catalog to the new method is a trivial task, it requires at least 30 days 
of computation time. Therefore, the author pseudo-randomly chose 1500 objects from the 
NORAD satellite catalog for analysis. This chapter can simply be divided into four dis¬ 
tinct parts. The first part, Section 3.1, discusses the data gathering, SGP4 propagation, and 
contains a quick overview of the whole process. The second part includes all sections from 
Section 3.2 to Section 3.6, which can be conveniently named as the KAM theory building 
part. The third part, starting in Section 3.7, introduces the evolutionary steps to the desired 
KAM torus by differential correcting the initial conditions and the frequencies of the newly 
built torus using SGP4 and TLE data. The fourth part is a summary. 

3.1 Data Gathering and Overview of the Method 

This effort uses the TLE obtained from the website www. space-track. org. This site 
allows users to download the TLE of most satellites. The previous work by Wiesel and oth¬ 
ers showed that there are some properties of satellites which cause difficulties in KAM tori 
construction process. The author also made an initial test in order to define an optimal ec- 
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centricity threshold level because this effort is intended for low eccentric orbits. Although 
the author defined an eccentricity of 10 3 and small as optimal eccentricity for the pro¬ 
gram to converge, he will further analyze the highest limit for the eccentricity to converge 
by modifying the code. It is also known that there are no nearly circular periodic orbits in 
the close vicinity of the critical inclination, which is 63.435°, or the complementary angle 
of 116.565° [58]. Although there exist periodic orbits for polar orbits, the close vicinity 
of polar orbits is avoided because of numerically unstable Legendre polynomial recursions 
for the gravitational potential calculations. The correction is left for future studies. An 
altitude of 300 km and more is chosen as a criterion because air drag shrinks the size of the 
torus down in phase space. It is also known that orbital maneuvers destroy the KAM torus, 
but they are not listed as a criterion because maneuvers are unpredictable. Table 8 shows 
known issues in the previous works and the selection criteria of test objects for avoiding 
these issues. 

Table 8. Test Case Selection Criteria 


Issues 

Test Case Selection Criteria 

Air Drag 

Objects which have altitute of 300 
kilometers and more 

Resonance 

Objects with periods that are not 
nearly an integer multiple of Earth’s 
rotational period 

Critical inclination and 
Polar Orbits 

Objects which have inclinations that 
aren’t in the close vicinity of the criti¬ 
cal and polar inclinations. 

(0° < i < 60° and 67° < i < 87° and 
93°</< 113° and 119° < i < 180°) 


There are 7,938 space objects which match the criteria in Table 8. The author wrote a 
C++ script which retrieves 7,938 TLEs for a period of 2 months ,which is between Jan 1, 
2013 and March 1, 2013, from www.space-track.org server. This script also prepares 
essential input files that are required by low eccentricity KAM torus construction program, 
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TORIS 1 



Figure 16. The Orbit of Toris 1 Satellite over 2 Months 

which can be conveniently named as theory building program. The solar activity is rela¬ 
tively low between the selected dates. The dates for solar minimum can be selected, but 
the author intended to analyse the recently launched satellites as well. Then, 1500 TLEs 
are pseudo-randomly selected, and propagated for an orbit for a period of 2 months by 
Vallado’s revised SGP4 code [21]. Figure 16 shows the orbit of TORIS 1 satellite over two 
months, which is simply a torus shape. TLEs are propagated for an orbit period because it 
is known that the position errors are smallest some point within the data interval and grows 
with respect to time going outwards into the future, or back into the past [67]. It is also 
desirable that JSpOC releases the TLE with an epoch time within data interval. JSpOC 
fits the orbits for 3 to 4 days for LEO satellites and a couple of weeks for higher altitude 
satellites. Therefore, there is a progression of orbits, and that progression includes more 
information than each individual orbit has. The idea behind propagating each TLE for a 2 
months orbit is to obtain more accuracy through combining smaller chunks of less accurate 
propagation data. Because SGP4 and TLE prediction data are considered observational 
data, data from the calculation of dynamics are needed, which will be combined to form 
an estimate. The low eccentricity KAM torus provides the dynamical data. The low ec¬ 
centricity KAM torus construction starts with the calculation of the periodic orbit, which 
is basically a boundary value problem. Once the periodic orbit is built, it is transformed 
to Fourier series [23]. The solution for the periodic orbit in the zonal potential problem 
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yields two basis frequencies, which are the Keplerian frequency and the nodal regression 
rate. The other missing frequency, which is the advance of the argument of perigee rate, 
can be obtained from the Floquet solution of the periodic orbit. The Floquet solution also 
yields two time linear terms due to adjacent periodic orbits. The solution includes only 
the zonal potential so far. However, perturbations, which are air drag, and sectoral and 
tesseral harmonics perturbations, are added to the Floquet solution. Therefore, the problem 
becomes forced linear system problem. The order of perturbations suddenly drops from 
10 -3 to 10~ 6 , which satisfies the main dictate of KAM theorem of small perturbations bet¬ 
ter. In addition, the momenta are inertial velocity components defined in a rotating frame 
of reference because of the set up of the dynamics. Therefore, simply substituting veloc¬ 
ity components with momenta is a valid approach, see Section 3.2. Indeed, defining the 
forcing function as a function of Q\ and Q 2 in a reference frame that rotates with Earth 
where sectoral and tesseral perturbations are stationary leads to the usual set up of the fun¬ 
damental equation for the KAM theorem, see Equation 86. After the low eccentricity KAM 
torus is built, it is fitted to the observational data acquired from SGP4 and TLE prediction. 
Non-linear least square algorithm is used to form an estimate. The modal variables and the 
frequencies of the periodic orbit are updated for each iteration, see Section 3.7. Although, 
the author represents the process piece by piece, each main program is connected to another 
by GNU/Linux bash files, which are batch files for MS-DOS environment. Therefore, it is 
a simple matter to switch different parameters, such as air drag, and sectoral and tesseral 
perturbations, on and off. Moreover, human intervention is not required during the pro¬ 
cess due to the bash files which orchestrate the whole orbit determination procedure. This 
method is a numerically constructed perturbation theory which uses periodic orbits under 
the influence of only the zonal potential as the starting point [69]. 
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3.2 Satellite Dynamics 


This effort uses different reference frames, such as rotating and inertial frames of refer¬ 
ence. The inertial reference frame is the ECI frame, and the rotating frames have rotation 
rates , (0 — Cl where the frame of reference is tied to the regression of the ascending node, 
or co = to® where it is tied to Greenwich. In such a rotating frame of reference the inertial 
velocity vector becomes [69]: 


v — r — 


x— coy 
y + cox > > 
z 


(96) 


where co appears in the components of the inertial velocity because of transport theorem. 
Then the satellite’s kinetic energy becomes [69]: 


T = \ {(x~COy ) 2 + (y + COx ) 2 + (z) 2 ) 


(97) 


The generalized momenta can be found from the canonical Hamilton equations, pi — 
dT/dqu as [69]: 


P x = X - coy 

p y — y+ cox (98) 

Pz = Z 
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where the generalized momenta are the inertial velocity components. Then, the Hamilto¬ 
nian can be calculated from H = £,■ ptq — T + V as [69]: 


H = \ ( pI + Py+P 2 z ) + ®(yPx - XPy) 

-jiid-) Vi**) <"> 

n=0 m=0 V ©/ 

x (C nm cos mX + S nrn sin mX ) 

where /i is the Earth’s gravitational parameter, R- , is the radius of Earth at the equator, and 
C nm , S nm are the geopotential coefficients. This effort uses NASA EGM-96 gravity model 
with an order and degree of 20. The relationship between rectangular coordinates and the 
radius, r, the latitude, 8 , the longitude, X, can be shown as [69]: 


r = \Jx 1 +y 2 +z 2 


sin 8 
tan X 


z 

\/x 2 +y 2 

y 

X 


( 100 ) 


where the rectangular coordinates are in a reference frame tied to the Greenwich, which 
is the Greenwich referenced ECR frame, see Section 2.4. The expansion of the geopo¬ 
tential in the zonal harmonics, which is m = 0, is used for constructing the periodic orbit. 
The sectoral and tesseral harmonics are added as perturbations to the periodic orbit. The 
Hamilton’s equations can be represented as [69]: 


X 



( 101 ) 


where X is the physical state vector, X T = {.r.y.z, p x . p y . p-}. and Z is a 2n x 2 n matrix, 
where n is the number of coordinates, that conveniently arranges the Hamilton’s equations 
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in vector form as [69]: 


f * 1 


y 


x=l 



Py 


0 I 
-7 0 


dH 

dqt 

dH 

dp. 


= < 


dH 

dpx 

dH 

dp y 

dH 

dp z > 

dH 

dx 

dH 

dy 

dH 
dz j 


( 102 ) 


The linearization of the Hamilton’s equations is needed because this effort is an estimation 
problem. The equations of variation can be represented as [69]: 


d 2 H 

i=Z 3X2 =Ax ’ 


(103) 


where x is represented as small changes to the physical state vector X. 

Air drag perturbations are also added to the the combination periodic orbit and the 
Floquet solution. The momenta are inertial velocity components. Therefore, air drag ac¬ 
celeration can be represented as [69]: 

_ ' CdApo P / 2 | „2 | „2 n 

Cl drag ~ ~ ^ —— — yPx+Py+Pz P 

' , P ° (104) 

= -B'^\/pi + pj + p;p 

Ho v 

where B* is another measure for the susceptibility of a satellite to air drag , m is the satellite 
mass, and p is the atmospheric density at perigee. [60]. For further information, see 
Equation 84. 
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Symmetric About Polar Axis 



Figure 17. The Earth’s Zonal Gravity Potential 


3.3 Periodic Orbits 

The periodic orbit and Floquet theory are used as starting point for calculating the 
perturbations because the previous KAM torus construction methods failed for low eccen¬ 
tricity orbits and the periodic orbit resides at the core of the torus. The modified Laskar 
algorithm, developed by Wiesel, is one of the KAM torus construction methods. It yields 
three basis frequencies for orbits by Fourier series spectral analysis. The low eccentricities 
are hard to be detected for methods like Laskar frequency analysis method, which requires 
high eccentricities that cause more dramatic change in the frequency spectrum of the orbit. 
Although this fact is a disadvantage for the previous KAM torus construction methods, it is 
favorable for this effort because the structure lies at the core of the torus is a periodic orbit. 

In 1960’s, the periodic orbits were studied extensively because they are more realistic 
than 2BP, which describes every orbit periodic. The Newtonian point mass plus the zonal 
harmonics terms lead to the periodic orbit because the zonal potential about polar axis is 
symmetric. Figure 17 represents the symmetry of the Earth’s zonal potential about polar 
axis. Nearly circular orbits are periodic orbits in the Earth’s zonal gravity potential for 
every inclination, except in the close vicinity of the critical inclination [58]. The periodic 
orbits are closely related to the frozen orbits, which have nearly stationary eccentricity 
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Figure 18. The Regression of Orbit Plane in the Zonal Gravity Potential 

and argument of perigee. The repeated ground tracks are of great importance for mission 
planning. The periodic orbits defined in a reference frame tied to the rotating Earth have 
repeated ground tracks. Therefore, the periodic orbits have been studied as frozen orbits 
as well [44, 43]. It has been shown that nearly circular orbits are unstable at the criti¬ 
cal inclination, and the periodic orbits in the vicinity of the critical inclination have high 
eccentricities [17]. The numerical proof of the existence of periodic orbits with high ec¬ 
centricities has been done by Broucke [8]. Moreover, Wiesel has developed a new relative 
motion solution for satellites using the low eccentricity periodic orbits, but their potential 
use for orbit determination hasn’t been explored until he uses the low eccentricity periodic 
orbits to construct the low eccentricity KAM torus [69]. 

In the zonal problem, there are periodic orbits. However, there are a few difficulties in 
building the periodic orbit. It is known that there are two non-zero frequencies, which are 
the Keplerian frequency and the regression of the node. These frequencies aren’t known. 
The periodic orbit precesses at a certain rate which can only be determined when it fin¬ 
ishes its motion for one period, and the orbits aren’t periodic in inertial space. Therefore, 
defining the periodic orbit in a reference frame that rotates with the node is required. Fig- 
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ure 18 represents the regression of orbit plane under the influence of zonal potential. The 
periodicity condition needs to be set up in such a way that the nodal regression rate of the 
orbit plane is not needed. Therefore, the initial conditions are defined with position vector, 
r = {a' o . 0.0}, and momentum vector, p = {x () , v () cos i G , v () sin i () }. The orbital period, T, and 
inclination, i, define the orbit, while the unknowns are Z T = {x 0 .x 0 . v 0 }. Because Earth 
doesn’t have North/South symmetry, x a is required as a component of the momentum vec¬ 
tor, p. The orbit becomes periodic if the 3 conditions after one revolution, T, are set up as 
[69]: 

f ^( T ) ) 


G = 


r(t) • p(t) -x 0 x 0 
r(T)-r(r )-x 2 0 


= 0 , 


(105) 


where the conditions dictate the satellite to cross the equatorial plane, z(t) =0, to have 
the same radial velocity as at the beginning, r(r) -p(r) — x 0 x 0 = 0, and to be at the same 
distance from Earth as at the beginning. The actual orbit precession isn’t known, but can 
be calculates as [69]: 


1 _i f(t) -r(0) 1 

-COS , 7 — , . , = -COS —— 

|r(T)|jr(0)| T *(0) 


(106) 


A linearization of Equation 105 using Taylor’s series expansion is required to obtain 
the periodic orbit [69]: 

G = G(S) + — 5S = 0 (107) 

do. 

The expansion using the state transition matrix <E(t, 0) , which is calculated from parallel 
numerical integration with the trajectory, yields the derivative matrix of G [69]: 


dG 

~dZ 


d G 


dX(z) 

dG 


dX( T) + 


dZ 


dX{ T 


-d>(T,0) + 


dG dX{ 0) 
dX{(y) dz 
dG \ dX{ 0) 


dX(0) J dz 


(108) 
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where <F(t, 0) can be calculated from <9 X(t)/< 9X(0). The partial derivatives in the Equa¬ 
tion 107 can be computed as [69]: 


dG 

<9X(t) 


= < 


0 0 1 

x y z 

2.x 2 y 2 z 


0 0 0 

x y z > 

0 0 0 


dG 

r)X(O) 


0 

-i 0 


-2x 0 


0 0 
0 0 
0 0 


0 0 0 
—x 0 0 0 > 
0 0 0 


dx(0) 

dZ 


= < 


10 0 0 
0 0 0 1 
0 0 0 0 


T 

0 0 

0 0 > 

cos i a sin i Q 


(109) 


Then, the initial condition parameter vector, E, can be corrected by the correction, <5E, 
which can be calculated by Equation 107 [69]. 

After the periodic orbit is built, it is represented as a Fourier series by harmonic analysis 
[23]. Figure 19 shows the representation of the periodic orbit by two angle variables. In 
Figure 19, Q\ is ’’mean argument of latitude analogue”, and (?2 is regression of the right 
ascension of ascending node. Q\ and Qi are the physical coordinates of the KAM torus 
[69]: 

£i=M + Q, (110) 


Q 2 = n~d g , ( 111 ) 

where M is mean anomaly, Q. is right ascension of ascending node, and 0„ is angle between 
the vernal equinox and the prime meridian. The periodic orbit is a function of only Q\ in a 
frame of reference, which is the ECNF frame, see Section 2.4, that rotates with the ascend- 
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Figure 19. The Representation of the Periodic Orbit by Angle Variables 


ing node. The nodal orbital frequency, C0i, and the inertial nodal regression rate, Oh, can 
be obtained from the Fourier series. These Fourier series store only the coordinates, , q, 
and £, because the momenta can be obtained by the kinematic relations in Equation 98. 
Moreover, it is a simple matter to generalize from ECEF frame to ECI or ECEF frames 
because the zonal potential is symmetric about polar axis. The inertial nodal regression 
rate, Oh, needs to be updated to define the regression rate with respect to the Greenwich. 
Equation 112 represents that the periodic orbit defined in ECNF frame, where the periodic 
orbit, X, is a function of only <2i, can be transformed to ECI or ECEF frame with a double 
rotation matrix about the z axis by the nodal angle <22 [69]: 




Rz(-Qi) 0 
0 Rz(-Qi) 


( 112 ) 


Equation 112 can also be used when the desired reference frame is ECI instead of ECEF 
with small modifications. Both <22 in the equation are substituted by the usual initial frame 
node, Q., and Oh is reduced by the Earth’s rotation to transform the periodic orbit, X, from 
ECNF frame to ECI frame. Although the two frequencies have been obtained from the 
periodic orbit, the third frequency is required for the solution. 
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3.4 Floquet Solution 


Periodic orbits plus the Floquet solution provides a method that numerically solves dy¬ 
namical problems which don’t have closed form solution. The linearization of the periodic 
orbit describes the behaviour of the system near the periodic orbit. Linearizing about the 
periodic orbit results in a linearization problem [64]: 

x = A{t)x , (113) 

where x is the first order displacements from the state X. Equation 113 is a time periodic 
linear differential equations. The solution to time periodic linear differential equations 
yields stability information of periodic orbits. The solution to the Equation 113 can be 
obtained as [64]: 

®(t,t 0 ) = E(t)exp(J(t -to))E~ l (t 0 ), (114) 

where E is a periodic matrix, J is a matrix of system frequencies in the Jordan normal form. 
The system frequencies are also called Poincare exponents. 

The Floquet solution can be obtained by calculating the Jordan normal form, J, and the 
periodic matrix, E, over one period. Over one period, Equation 114 becomes [64]: 

<J>(t,0) =E(t)exp(Jt)E~ l (0), (115) 

where E(t) = E( 0) because E is periodic. Figure 20 represents the visual depiction of 
E periodic matrix. Geometrically, defining E{t) — E(0) puts a sets of coordinates on the 
periodic orbit, and they follow the trajectory by rotating around themselves. At the end 
these sets of coordinates smoothly join themselves. is the state transition matrix over a 
period. Indeed, because the system is periodic, we have <4> for all time. The <t> matrix is 
called monodromy matrix. The eigenvalues and eigenvectors of the <t> matrix can be found 
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E( 0) = E(t) x 


Figure 20. E Periodic Matrix 


by rearranging Equation 114 [64]: 




(116) 


where the exponential of the Jordan form ,J, is the eigenvalue matrix of <t>. E periodic 
matrix yields the eigenvalues of the <t> matrix. If A/ are eigenvalues of the monodromy 
matrix, Equation 116 becomes [64]: 



(117) 


where ft), are the system frequencies, which are also called Poincare exponents. The inter¬ 
pretation of the Poincare exponents is similar to the eigenvalues of a constant coefficient 
system. For Hamiltonian systems without dissipation, the frequencies must be a pair of 
positive / negative real numbers or a pair of a pair of imaginary numbers. As an exception, 
when A, = 0, there is a repeated value of ft)/ = 0. Repeated roots in linear systems lead to 
off-diagonal time terms. These pair of zeros are called a degenerate mode. For this effort, 
the Hamiltonian has two integrals of the motion, which are the Hamiltonian itself, and the 
z component of the angular momentum. Because each of these pair of zeros correspond 
to one integral of motion, this effort has two degenerate modes. Each degenerate mode 
involves a static displacement and a linear drift is associate momentum is not zero. These 
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static displacements can occur either along the orbit or around the polar axis. These two 
degenerate modes correspond to two frequencies. The third frequency is the free eccen¬ 
tricity of the orbit, which is a pair of imaginary frequencies by the eigenvalue calculation, 
because the imaginary part of <y, is the oscillatory frequency of the mode i. For this effort, 
the oscillatory frequency is ±(03 with algebraic sign that depends on the inclination. As a 
result, the calculation of eigenvalues yields the complete linear stability information of the 
periodic orbit [64, 69]. 

The calculation of eigenvector yields E(t) over one period, which can be used to trans¬ 
form from physical variables to model variables. The eigenvalue problem for equation is 
highly singular, and a special treatment is required. The differentiation of the state transi¬ 
tion matrix becomes [63]: 


Mt,to) =E(t)e J{t - t ° ) E~\t 0 )+E{t)Je J(t - t ° ) E-\t 0 ) (118) 


By substituting into Equation 119 : 

<P(t,t 0 )=A(t)<P(t,t 0 ) (119) 

A differential equation for F(t) becomes: 

E(t)e J ^- to) E-\t 0 )+E(t)Je J ^- to) E-\t 0 ) = A[t)E{t)exp{J[t - t 0 ))E~ l (t 0 ) 

E(t) +E(t)J —A(t)E(t) (120) 

E —AE — EJ 

Since E( 0) is obtained from the periodic orbit, solving Equation 120 from / = 0 to / yields 
F(t) over one period. E(t ) periodic matrix is reduced to Fourier series [23]. For this effort, 
the periodic matrix is a function of Q\, E(Q\). 

The transformation from physical coordinates to modal variables is possible because 
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F(t) periodic matrix is calculated. Assuming y is the modal variable [63]: 


y = E~\t)x 


(121) 


By substituting into Equation 122 : 


x(t) = <j>(t,t 0 )x(t 0 ) 


( 122 ) 


The Floquet solution represented in modal variables becomes : 

Table 9. Analogous of Modal Variables 


Modal 

Variables 

Analogous of Modal Variables 

yi 

a displacement along the argument of latitude di¬ 
rection, Q\. 

yi 

a change in the momentum which is analogous to 
the Delaunay momentum, L. 

T3 

an offset in the nodal variable, Qi- 


a change in the z component of orbital angular mo¬ 
mentum. 

T5 

ecos co 

ye 

esin co 


x{t) = E{t)e J ^E-\t 0 )x{t 0 ) 

E~ l {t)x{t) = e J ^~ t °^E~ 1 (t 0 )x(t 0 ) (123) 

y (t) = e J{t - to) y(t 0 ) 


Equation 123 is a solution to the linear system to the form [69]: 


y = Jy 


(124) 
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The Floquet solution can be added to the periodic orbit solution as : 

Xecr = Rz(-Qi) {Xpo(Qi) + E(Qi)exp(J(t - t 0 ))y(t 0 )}, (125) 

where double z rotation is represented as R z , the nodal frame periodic orbit series are shown 

as ( X)po . The modal variables , y\ and V 2 , are local versions of the global angles , Q\ and 

Q 2 , and they are scaled to suit this purpose. The Jordan form matrix becomes : 

' ' 

It —0 0 0 0 

0 1 0 0 0 0 

0 0 1 t — t 0 0 0 

exp(J(t — 1 0 )) — (126) 

0 0 0 1 0 0 

0 0 0 0 cos <23 sin Q 3 

0 0 0 0 —sin Q 3 cos Q 3 

The modal variables are analogous to 6 orbital elements that can be used to represent 
nearby motion. Table 9 represents an analogy between modal variables and miscellaneous 
orbital elements. The global variables, Q\ and Q 2 , should be freely traded with the modal 
variables, y\ and y 3 , because this will keep them from growing unboundedly [69]. 

3.5 First-Order Perturbations 

The periodic orbit and its Floquet solution has been obtained from Equation 125. This 
solution has the Newtonian term and all of the zonal gravity harmonics. The perturbations 
start approximately on the order of 10 5 instead of 10 3 , which is advantageous because 
errors originated from the calculation of perturbations are less influential. However, the 
method for calculating first-order perturbations can be extended to higher order perturba¬ 
tions. This effort includes only sectoral and tesseral harmonics and air drag. However, the 
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exact same method applies to lunar and solar mass perturbations by including at least one 
extra angle in the Fourier series which describes the motion of the Moon or the Sun with 
respect to Earth in a frame that rotates with it. The perturbing force, X pert , can be added to 
Equation 113 as [69]: 

X = Ax + Xpert (127) 


The perturbing acceleration can be expanded about the periodic orbit because this is a 
forced linear system [69]: 
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(128) 


The first-order perturbations can be added to the Floquet solution by Equation 124 as 
[69]: 
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0 0 0 0 - ft> 3 0 


(129) 


where for sectoral and tesseral perturbations, the trailing term, E~ l X pert , is a function of 
the global variables, Q\ and Qi- The forcing term can be transformed to a double Fourier 
series in these angles. This applies to air drag as well. However, the air drag forcing term 
is a function of only one global variable, Q\. 

For the first-order perturbations, a perturbing acceleration transformed to the modal 
variables is decoupled. Two different terms, which are the constant and periodic Fourier 
series terms, are examined for the perturbing acceleration as the response to the oscilla- 
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tory and the two degenerate modes. Each factor is evaluated separately. A constant solu¬ 
tion is obtained from the constant terms in the Fourier series for the oscillatory mode by 
Equation 129 as [69]: 
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Then, the solution to the Equation 130 becomes [69]: 
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(131) 


These terms can’t be ignored because (O 3 is a small number for the most part. The modal 
variables, y\ and y 2 , don’t become infinite near the periodic orbit, which has zero eccen¬ 
tricity for this effort, because they are analogous to the quantities, ecos co and esin co, in 
the 2BP. The constant terms in the Fourier series for one of the degenerate modes doesn’t 
yield a constant solution [69]: 




Then, the solution becomes: 
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( yi \ ( (ci-y 2 (t. 0 )){t.-t 0 ) + \c 2 (t 2 -t^)\ 


w 
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(132) 


(133) 


where the solution for y 2 = c 2 is written at the start time because the other secular terms, 
which appear due to the periodic perturbations in y 2 , can be absorbed. The solution for y\ 
has been obtained by substituting >’2 into yi. Equation 133 represents that y 1 and V 2 could 
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have quadratic terms, and y 2 and 34 could have a linear rate of change. When there is no 
dissipative force in the system, C2 should be zero, and a quadratic term isn’t expected for in¬ 
track and nodal perturbation. For air drag perturbation, quadratic behavior is expected.)’1 
and y 2 can be observed by their associated global variables because the modal variables 
aren’t scaled. The exact same approach is taken for periodic Fourier terms. 

The periodic terms in the Fourier series for a degenerate mode yields the forced linear 
system as [ 69 ]: 
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C 1 cos{niQi+n 2 Q2)+si sin(niQi + n 2 Q 2 ) 
c 2 cos(niQi+n 2 Q2) + s 2 sin(n\Q\ + n 2 Q 2 ) 


\ 


( 134 ) 


where the setup is for sectoral and tesseral harmonics, but if a perturbation requires it, 
an additional angle can be introduced. In the equation, n\ and n 2 are integers, and w p = 
n\( 0 \ -\-n2O1h. The solution for 34 can be calculated [ 69 ]: 


yi(t) = (— sin(n\Q\ +112Q2) 
(Op 

Then, substituting into the y\ yields : 


— cos(niQi +n 2 Q 2 )) 
(Op 


( 135 ) 


y\(t) = <( (— - %) sin(niQi -\-n 2 Qi) ~ (— - %) cos(n\Q\ +n 2 Q2) 

'On cot: COr, (Q 


Q2 S 2 

- (— sin(niQi +n 2 Q2) -- cos(rnQi -\-n2Q2)) 

(On (On 


( 136 ) 


where the constant terms in the y 2 solution lead to time dependent terms in the solution 
for yi. However, the secular terms introduced by the constant terms in the y 2 solution has 
already been absorbed into the solution provided by Equation 133 . 
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The remaining forced linear system is due to the periodic terms in the oscillatory mode. 


Equation 137 represents the forced linear system [69]: 
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Let both a forced solution and the perturbation term has the same form: 
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(138) 


After inserting assumed solution and simplifying it, the solution becomes a set of linear 
equations for the forcing coefficients as: 
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(139) 


Once the perturbing Fourier series is obtained, it is converted to the perturbation solution, 
and stored. This method can be applied to any other perturbations with simple modifica¬ 
tions, and the code written for this method requires only the new data file and truncation 
level associated with the perturbation. It should be noted that Equation 139 includes the 
expected small divisors problem, which will lead to problems near the resonance. 

Equation 125 can be modified to include the first-order perturbations as[69]: 


Xecr —Rz{~Qi) {Xpo(<2i) +E(Q\) [exp(J(t — t 0 ))yf ree ifo ) + }'forced\ } > (140) 
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where the modal variables from the Floquet solution is represented by yf ree , and the per¬ 
turbed solution is denoted by yf orcee j . The perturbed solution, yforced’ can be a function of 
different angles, depending on the perturbation. 


3.6 Second-Order Eccentricity Perturbations 


The modal variables, yi and y$, are necessary for fitting orbit and for small perturba¬ 
tions. The global angles Q\ and Qi absorb the large changes. Moreover, secular terms may 
be absorbed in their frequencies, (0\ and co 2 . However, the expansion of the torus along, V 5 
and y 6 , which describe the eccentricity and argument of perigee states, is necessary because 
not all orbits have small enough eccentricities. The equation for the expansion of the two 
model variables, y$ and yg, when setting the other model variables approximately zero in 
physical variables becomes [69]: 


Xi - Zr'< 


d 2 H 


dXadXp 


Xpo 
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d 3 H 

X P + 2 Zia dX a dXpdXy 


Xfi Xy+ ... 


Xpo 


= Aia(Q\)x a + xBjafl (Ql)xccXp 


(141) 


where Z is the symplectic matrix, i = (1...6) for this specific case, the linear, time varying 
differential equations are the equations of variation. Representing the equation in the modal 
variables, Equation 141 becomes: 


1 

yi = Jiayoc + xE ia B a p y Ep 8 E ye y 8 y e + ... 

\ (142) 

= Jia,ya + 2^1apyuyp + ••• 

where 6 '^, the periodic orbit, and the E matrix are periodic functions of Q\. The quadratic 
matrix, B\- k with y^ and yg portion in the final two indices yields a periodic six by two by 
two tensor [69]. 

Assuming all model variables are zero except the eccentricity mode, the linear matrix 
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becomes: 



Then, the linear solution can easily be obtained as: 


(143) 


y{t) =®(t,t 0 )y{t 0 ) = 


cos(0"i (t — t 0 ) sinoy$(t — t 0 ) 
—sinGh,[t — t 0 ) cos(Q?,(t — t 0 ) 



(144) 


If the quadratic terms in the Equation 142 is small enough, the first order solution, y(t) — 
<&(t,t 0 )y(t 0 ), can be used for the equation as: 


>’( — ha y a + fiSyyityy s{to) 

f (145) 

= hay a + 2 B,, iap y a(t 0 )yp (t a ) 

where B ” is a function of both Q\ and Q 3 = CO 3 (f — t 0 ). The solution to Equation 145 is 
assumed to be generalized as [69]: 


yiit) = d>i o 1] y a (t 0 ) + ^^pyaMypito) + ... , (146) 

where d> (1 )(^ 0 ) = /,as it is in the classical linear system theory, and d>( 2 )(t 0 ,t 0 ) = 0 be¬ 
cause the initial conditions, yi(t 0 ), need to be preserved. Therefore, the time derivative 
of Equation 146 can easily be calculated because yi(t 0 ) are constant. Inserting the time 
derivative in the form of the modal differential Equation 145 yields [69]: 

y i = M! ] ya(t 0 ) + ^apyaMypito) + - 
= J ifMp a ya{ta) + ^JiM/apy^yp (to) ( 147 ) 

+ ^ B ”iapya(to)yp(to) 
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where the arbitrary indices are for a uniformity requirement in the initial conditions. The 
first order state transition matrix differential equation can be obtained as [69]: 


(148) 

Then, the new result in the second order becomes: 

7*. (M9) 

where it is already known that < $> i j k {t 0 .t 0 ) — 0. The solution for Equation 149 has been 
provided for the constant terms in the Fourier series of B” by Equation 131 and by Equa¬ 
tion 139 for the periodic terms. If indices j and k are constant, the solution to the system 
will be same as the first order case. Therefore, O' 21 will be a periodic function of the global 
variables, Q\ and Q^. The extension to higher order of the eccentricity/argument of perigee 
is believed to better characterize the short periodic oscillations. 

3.7 Least Squares Fit to SGP4 and TLE Data 

For this method, the point on the periodic orbit that corresponds to the epoch time of 
SGP4 and TLE data should be calculated in terms of Q\ and Qi- From the perspective of 
estimation theory, a reference trajectory should be corrected to yield an estimate of the ac¬ 
tual trajectory. For this effort, the low eccentricity KAM torus predictions are considered as 
the reference trajectory, and the SGP4 and TLE predictions are considered as the observa¬ 
tional data. Wiesel and Frey showed that the KAM torus basis frequencies can be obtained 
from SGP4 and TLE data. The desired periodic orbit is the one with zero y 2 and 34 model 
variables. The least squares routine for this effort fits the low KAM torus to the SGP4 and 
TLE data by updating the global variables, Q\ and Q 2 , the frequencies, 0 )\ and ah, and B*. 
The convergence criterion is met when the estimate correction vectors of modal variables, 
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Error Ellipsoid 

Figure 21. The Converged and Unconverged State Correction Vectors 

8yt, are smaller than 1% of the sigma values, ay Convergence for least squares method 
depends on the correction vectors and the error ellipsoid. If the correction vectors lie well 
within the 1 a error ellipsoid [67]. Figure 21 represents the converged and unconverged 
state correction vectors. 

This effort uses nonlinear least squares. This method combines the SGP4 and TLE 
data, which has been propagated for an orbit over two months, and the low eccentricity 
KAM torus to form an estimate. For each observation, 6, a series of calculations are made 
for nonlinear least squares method. First, the state transition matrix is obtained , <t>(tj,t 0 ), 
by propagating the state vector to the observation time, tj. Second, the residual vector, 
r,- = Zi — G(x), is calculated. The SGP4 and TLE data are represented as zu and the low 
KAM torus data are represented as G(jc). Third, the observation matrix, T, = //,<$>, the 
covariance matrix, P§ x , and the state correction vector at epoch are calculated as [67]: 

Tj = Hj<£> 

Ps, = CLPQT't ,)- 1 (150) 

i 

Sx(to ) = Pdx^T? Qr'n 

i 

where Q\ is a matrix that has sigma errors associated with the measurements, H l is obser¬ 
vation function linearization matrix. Then, the reference trajectory is corrected as: 

%ref+l do) — Xrefifo) T 8x(to) (151) 
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Table 10. Three Shortcomings of the Least Squares Method [11] 


1 

Although every observation has different accuracy for the 
observations, each observation is weighted equally, which is 
a problem when there is observation with poor accuracy be¬ 
cause the least squares method will be biased towards them. 

2 

The errors can be correlated, not independent as least 
squares method considers them to be. 

3 

The method soesn’t consider that the errors are samples 
from a random process, and it doesn’t utilize any statisti¬ 
cal information 


A 



NULL SPACE OF MATRIX A COLUMN SPACE OF MATRIX A 

Figure 22. The Visual Depiction of the Least Squares Method 


Finally, the process ends when a desired limit for the convergence is met. However, there 
are short comings of the least squares solution, which are listed in Table 10. In addition, 
Figure 22 shows the visual depiction of least squares method. The least squares method 
yields the the projection of the solution in the column space, where the solution is supposed 
to exist. Therefore, in a sense, the least square method minimizes (Ax — b ) 2 , and yields the 
closest approximate of the solution. 


3.8 Summary 

The geometry of the core of torus in the full geopotential appears to be a multiply- 
periodic two dimensional ribbon structure as a function of the argument of latitude angle, 
<2i, and the nodal angle, Qi, whereas it is a periodic orbit in the zonal potential. The orbital 
eccentricity is a mode of oscillation, which has been expanded to the second order. For 
the low eccentricity KAM torus theory, each satellite has a set of data files for the periodic 
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orbit, the modal E matrix, and perturbation series in the form of Fourier series. The periodic 
orbit has two large angles, Q\ and 02. and four small parameters, >’2, V4, V5. and y$. The 
change in the two angles is quite large, but the small parameters are expected to be small all 
the time. The theory includes two parts. The first part is a data package of Fourier series, 
which need occasional update. The second part is a set of seven parameters, Q\, Q 2 , y 2 , J 4 , 
y 5 , yg, and B*, at epoch t 0 which require more frequent update. Multiply-structured Fourier 
series define the trajectory under modeled perturbations, and it is specialized to a particular 
satellite because it is a numerical method. However, it takes less than a minute on a modern 
personal home computer [69]. After the theory files are built and the SGP4 and TLE 
propagation is obtained, the low KAM torus is fitted to the SGP4 and TLE data because 
the point on the torus that matches the epoch time of SGP4 and TLE data must be found. 
Moreover, the fitting process corrects the frequencies of the torus, (0i and ah, and it figures 
out the B * value for the propagation period. It should be noted that the actual accuracy 
that can be achieved by this theory is supposed to be far better with the raw observational 
data than SGP4 and TLE data due to the inaccuracies in the SGP4 and TLE. The whole 
process from the SGP4 and TLE propagation to the least squares fitting needs no human 
intervention because all scripts are linked by bash files. The author modified every script 
to be both compatible for Windows™, and GNU/Linux platforms. The computational part 
of this effort has been conducted on a 64-bit GNU/Linux machine. 
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IV. Results and Analysis 


In this chapter, the results will be presented and discussed. 1500 TLEs were obtained 
from www. space-track. com servers, and each were propagated for an orbit for 2 months. 
Then, the low KAM torus theory files were created, which are associated with the periodic 
orbit, sectoral and tesseral, eccentricity, and air drag perturbations. Finally, the low eccen¬ 
tricity KAM torus was fitted to the SGP4 and TLE propagation data, and the residuals were 
calculated. Therefore, the least squares fitting routine not only locates the epoch point on 
the KAM torus and corrects the frequencies, but also yields residuals. The first section 
presents the success rate for convergence between the low eccentricity KAM torus and the 
SGP4 and TLE data, and provides details related to the root causes of failures. Each failure 
was analyzed individually. The second section presents the relationship between orbital 
characteristics and residuals for all 1500 test cases. Again, some plots for failed cases are 
provided to figure out the root causes of the failures. The third section discusses and ana¬ 
lyzes position and velocity residual plots for LEO and GEO objects in the RAN reference 
frame, see Section 2.4. The fourth section describes the optimal mean orbital limitations 
that never fail the convergence of the low eccentricity KAM torus. The fifth section com¬ 
pares the low eccentricity KAM torus and the SGP4 predictions for the best fit and a mean 
test case. The last section summarizes the chapter. 

4.1 Overview of Results 

All publicly available TLEs which have eccentricity on the order of 10 3 and smaller 
were retrieved from the www. spacetrack. com servers. 1500 TLEs were pseudo-randomly 
selected to fit the conditions mentioned in Section 3.1. First, each TLE was propagated for 
an orbit for 2 months period of time. Next, the periodic orbit, and its linearization, which 
is Floquet solution, for each test case was calculated, and perturbations were added to the 
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Table 11. Rms Values for Least Squares Fitting of 1500 Test Cases 


Rms range (Km) 

Percentage 

0-0.5 

0.6% 

0.5-1 

8.53333% 

1-3 

8.13333% 

3-5 

3.73333% 

5-10 

5.53333% 

10-20 

8.93333% 

20-40 

8.06667% 

40-100 

4.86667% 

100-200 

4.2% 

>200 

4.26667% 


solution. Then, the SGP4 and TLE data were fit using least squares, and the residuals for 
each test case were calculated. The success rate for convergence is 56.8667%. Table 11 
represents the percentage of different range of rms values for all 1500 test cases. It should 
be noted that these percentages will increase if the propagation time period is smaller, be¬ 
cause air drag is stochastic, and 2 months is long enough to terminate the torus. Therefore, 
another analysis was conducted. One hundred failed test cases were pseudo-randomly se¬ 
lected and propagated for 20 days instead of 2 months. The success rate of least squares 
fitting over 20 days is 66%. Table 12 shows the rms values for the 100 previously failed test 
cases. It is clear that the air drag is the root cause of the failure for 66 test cases. The other 
34 test cases will be analyzed individually in order to find the root causes of the failure. 
Table 13 represents the 34 failed test cases for the second time. 

Iridium 68, NORAD id 25291, is an operational satellite. Therefore, maneuvers for 
orbit maintenance, or attitude control maneuvers might cause the failure. Another reason 
for the failure might be its relatively low height of semi-major axis and high B* values, 
which are 7.1558000c+ 03 km and 2.0938000c — 04, respectively. Iridium 68 was propa¬ 
gated from 02-01-2013 to 02-09-2013 in order to check whether the air drag effect is the 
reason for the failure. However, the fitting process yields residuals on the order of 10 7 for 
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Table 12. Rms Values for Least Squares Fitting of the 100 Previously Failed Test Cases 


Rms range (Km) 

Percentage 

0-0.5 

0% 

0.5-1 

3% 

1-3 

15% 

3-5 

5% 

5-10 

9% 

10-20 

3% 

20-40 

5% 

40-100 

4% 

100-200 

8% 

>200 

14% 


the second time, which is a typical error characteristic for polar inclination. Therefore, it 
is concluded that its inclination of 8.6397400e + 01° causes the failure. In addition, all the 
Iridium satellites are in a Walker constellation, and they have nearly the same inclinations. 
However, the author tested all other Iridium debris and satellites to confirm whether the 
high inclination is the reason for the failure. 

Iridium 33, NORAD id 35078, is a debris. Its inclination, semi-major axis, and B* 
values are very close to what Iridium 68, NORAD id 25291, has. Iridium 33, NORAD id 
35078, was propagated from 03-01-2013 to 03-08-2013 in order to check whether the air 
drag effect is the reason for the failure. However, the residuals are identical to the residuals 
that Iridium 68 yielded during the fitting process. Therefore, it has been concluded that 
its inclination of 8.6399100c+ 01° causes the failure. Moreover, there are 4 other Iridium 
33 debris which have approximately the same inclinations. Thus, all Iridium 33 debris in 
Table 13 fail due to their high inclinations. 

Meteor 1-21, NORAD id 7714, is an operational weather satellite. Its period is 6.1373891^-1- 
03 seconds, which is in the close vicinity of the 14:1 resonance. It was propagated from 09- 
02-2013 to 09-08-2013. It yielded residuals on the order of hundreds of meters during the 
fitting process although the convergence was accomplished. Because the fit was converged, 
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Table 13. 34 Test Cases Failed for the Second Time 


NORAD id 

Satellite Name 

Type 

Launch Date 

25291 

IRIDIUM 68 

PAYLOAD 

1998-04-07 

35078 

IRIDIUM 33 

DEBRIS 

1997-09-14 

7714 

METEOR 1-21 

PAYLOAD 

1975-04-01 

29349 

KOREASAT 5 

PAYLOAD 

2006-08-22 

25531 

IRIDIUM 83 

PAYLOAD 

1998-11-06 

26599 

BEIDOU 1 

PAYLOAD 

2000-10-30 

34356 

IRIDIUM 33 

DEBRIS 

1997-09-14 

16495 

COSMOS 1726 

PAYLOAD 

1986-01-17 

33697 

FENGYUN 1C 

DEBRIS 

1999-05-10 

9701 

DELTA 1 

DEBRIS 

1973-11-06 

21153 

SL-8 R/B 

ROCKET BODY 

1991-03-12 

24115 

PEGASUS 

DEBRIS 

1994-05-19 

23092 

COSMOS 2279 

PAYLOAD 

1994-04-26 

31039 

FENGYUN 1C 

DEBRIS 

1999-05-10 

15308 

GALAXY 3 

PAYLOAD 

1984-09-21 

28923 

FREGAT R/B 

ROCKET BODY 

1991-03-12 

31988 

FENGYUN 1C 

DEBRIS 

1999-05-10 

34928 

IRIDIUM 33 

DEBRIS 

1997-09-14 

30760 

FENGYUN 1C 

DEBRIS 

1999-05-10 

31561 

FENGYUN1C 

DEBRIS 

1999-05-10 

29830 

FENGYUN 1C 

DEBRIS 

1999-05-10 

31339 

FENGYUN1C 

DEBRIS 

1999-05-10 

33605 

FENGYUN 1C 

DEBRIS 

1999-05-10 

29921 

FENGYUN1C 

DEBRIS 

1999-05-10 

20061 

NAVSTAR 14 

PAYLOAD 

1989-06-10 

7560 

SL-8 DEB 

DEBRIS 

1972-08-16 

36482 

IRIDIUM 33 

DEBRIS 

1997-09-14 

13073 

COSMOS 1275 

DEBRIS 

1981-06-04 

24435 

EXPRESS 2 

PAYLOAD 

1996-09-26 

21782 

COSMOS 2168 

PAYLOAD 

1991-11-12 

14966 

SL-8 R/B 

ROCKET BODY 

1984-05-11 

35981 

COSMOS 2251 

DEBRIS 

1993-06-16 

34488 

IRIDIUM 33 

DEBRIS 

1997-09-14 

23404 

COSMOS 2297 

PAYLOAD 

1994-11-24 


86 




it was propagated for some random dates over 5 days in order to observe an improvement 
in the residuals. However, the residuals appeared to be poor for all attempts. Its semi-major 
axis and B* values can’t cause such poor fits. Therefore, it is concluded that the resonance 
is the reason for the poor residuals. The KAM torus might deform rapidly near the resonant 
orbits if not fail completely. 

Koreasat 5, NORAD id 29349, is an operational satellite. It was propagated from 05- 
01-2013 to 05-10-2013 and from 07-01-2013 to 07-08-2013 in order to check whether the 
third body mass perturbations lead to the failure. Both propagations yielded the same big 
residuals during the fitting process. Its period is 8.6165984c 4-04 seconds, and it is in the 
close vicinity of the 1:1 resonance. Therefore, it is concluded that the resonance is the 
reason for the failure. 

Iridium 83, NORAD id 25531, is an operational satellite. It was propagated from 02-01- 
2013 to 02-08-2013. The fitting process yielded characteristic polar inclination residuals, 
which were on the order of 10 7 km. Therefore, it is concluded that high inclination, which 
is 8.6396500c + 01°, is the reason for the failure. 

Beidou 1, NORAD id 26599, is an operational satellite. It was propagated from 01-01- 
2013 to 01-10-2013 and from 09-01-2013 to 09-10-2013. Both least squares fits yielded 
big residuals, and the convergence failed for both attempts. The resonance might be the 
reason for the failure because its period is 8.7352213c + 04 seconds. 

Cosmos 1726, NORAD id 16495, is an operational satellite. It was propagated from 
01-01-2014 to 01-14-2014. It has an altitude of 527 km, and aB* value of 3.4805000c —04. 
The orbit was successfully fitted to SGP4 and TLE with an rms value of 7.54832 km over 
13 days. The poor fit shows that the atmospheric model used in the program, which is 
based on the U.S. Standard Atmosphere, should be improved to predict satellites with low 
altitudes with better accuracy. Figure 23 represents the plots for position and velocity 
residuals. In the plot, along-track residuals are relatively bigger than the other components, 
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Position Residuals 



Time (Julian Dates) 
Velocity Residuals 



Time (Julian Dates) 


Figure 23. The Residuals for Cosmos 1726 over 13 days 

which proves that the air drag effect is the reason for the failure of the first attempt. 

Fengyun 1C debris, NORAD id 33697, 31988, 30760, 31561, 31339, 33605, and 
29921, are in the close vicinity of the 14:1 resonance. Each debris was propagated for 
a short period of time in order to figure out whether the air drag effect is the reason for the 
failure. None of the fits converged. Therefore, the reason for the failure is the 14:1 reso¬ 
nance. Moreover, Fengyun 1C debris, NORAD id 31039, and 29830, prove this statement 
because although two of them have lower height and higher B* value, both fits converged, 
and the residuals weren’t bad for their low altitude. Fengyun 1C debris, NORAD id 29830, 
yields an rms value of 9.47329 km, and NORAD id 31039 yields an rms value of 22.8057 
km. Figure 24 represents the plot for position and velocity residuals of Fengyun 1C de¬ 
bris, NORAD id 29830, and Figure 25 shows the plot for position and velocity residuals of 
Fengyun 1C debris, NORAD id 31039. In the plots, the quadratic characteristic of air drag 
effect can be easily seen. 

Delta 1 debris, NORAD id 9701, has relatively high B* value, which is 1.2101000c — 















Residuals (Km/s) 


Position Residuals 


10 

5 


- 5 

!2 -10 
Q: -15 
-20 
-25 

23012 23013 23014 23015 23016 23017 23018 23019 

Time (Julian Dates) 

Velocity Residuals 

0.015 
0.01 
0.005 
0 

-0.005 
-0.01 
-0.015 
-0.02 
-0.025 
-0.03 
-0.035 

23012 23013 23014 23015 23016 23017 23018 23019 

Time (Julian Dates) 


* 

b 

1 

* 1 1 

Radial (U) * 

Along-Track (V) ■ 

Cross-Track (W) ■ 

1 





Figure 24. The Residuals for Fengyun 1C, NORAD id 29830, over 6 Days 
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Figure 25. The Residuals for Fengyun 1C, NORAD id 31039, over 6 Days 
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02, and it is also in the close vicinity of the 13:1 resonance. It was propagated from 01- 
01-2014 to 01-09-2014. The least squares fit converged with an rms value of 94 km over 8 
days.Therefore, the air drag effect is believed to be the prime reason for the failure. 

SL-8 R/B, NORAD id 21153, is a rocket body. It has an inclination of 8.2926500c + 
01°, a semi-major axis of 7.3573300e + 03 km, a B* value of 1.3800000c —04, and a period 
of 6.2804555c+ 03 seconds. It was propagated from 01-02-2014 to 01-08-2014. The least 
squares fitting process converged with an rms value of 47.0068 km. Thus, the air drag is 
the reason for the failure. Figure 26 represents the plot for position and velocity residuals. 
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Figure 26. The Residuals for SL-8 R/B, NORAD id 21153, over 10 Days 

Pegasus, NORAD id 24115, is a debris. Its relatively low altitude, which is 512.275 
km, and relatively high B* value, which is 6.3281000c — 04, are responsible for the failure 
because the least squares fit converged when it was propagated from 01-02-2014 to 01-09- 
2014 over 8 days. Figure 27 represents the plot for position and velocity residuals. 

Cosmos 2279, NORAD id 23092, is an operational satellite. It was propagated from 
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Figure 27. The Residuals for Pegasus, NORAD id 24115, over 9 Days 

01-01-2014 to 01-19-2014. The least squares fit converged with an rms value of 13.8748 
km over 18 days. Therefore, it is clear that the least squares fit has failed to converge due 
to the air drag effect. Figure 28 represents the plot for position and velocity residuals. 

Galaxy 3, NORAD id 15308, is an operational satellite with a period of 24.03. The 
propagation from 01-01-2014 to 01-10-2014 yielded an rms value of 1.94622 km. There¬ 
fore, the inaccuracies in the TLEs caused the failure. 

Fregat R/B, NORAD id 28923, is a rocket body. It has an inclination of 5.6419100e + 
01°, a semi-major axis of 2.9813700e + 04 km, and a period of 5.1231209e + 04 seconds. 
It was propagated for different time periods, and all attempts yielded characteristic rms 
values for resonant orbits, which were approximately on the order of 10 6 , and the residuals 
were random-like. Therefore, it is concluded that the resonance is the reason for the failure. 

Navstar 14, NORAD id 20061, is an operational satellite. It has an inclination of 
5.5107800e + 01 °, a semi-major axis of 2.7647OOOe + 04 km, and a period of 4.5749126e + 
04 seconds. It was propagated from 01-01-2014 to 01-10-2014. The least squares fitting 
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Figure 28. The Residuals for Cosmos 2279, NORAD id 23092, over 18 Days 


process yielded random big rms values for iterations, which were on the order of 10 6 . 
Therefore, it is concluded that the resonance is the reason for the failure. 

SL-8 debris, NORAD id 7560, has a period of 6.1487649c+ 03 seconds. Thus, it is in 
the close vicinity of the 14:1 resonance. It was propagated from 01-01-2014 to 01-10-2014. 
The least squares fit failed to converge, and yielded random big rms values for iterations, 
which is typical residuals for orbits in the close vicinity of the resonance. 

Cosmos 1275, NORAD id 13073, is a debris. It was propagated from 01-02-2014 to 01- 
09-2014. The convergence was accomplished with an rms value of 3.39862 km. Therefore, 
the reason for the failure is the air drag effect. Figure 29 represents the plots for position 
and velocity residuals over 7 days. 

Express 2, NORAD id 24435, is an operational satellite. It has an inclination of 
1.2856800c + 01°, a semi-major axis of 4.2182300c+ 04 km, and a period of 8.6219633c + 
04 seconds. It is in the close vicinity of the 1:1 resonance. Express 2 was propagated from 
01-03-2014 to 01-08-2014 over 5 days. It yielded random big rms values during the least 
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Figure 29. The Residuals for Cosmos 1275, NORAD id 13073, over 7 Days 

squares fitting process, and it failed to converge. Therefore, the resonance is the reason for 
the failure not the third body mass perturbations. 

Cosmos 2168, NORAD id 21782, is an operational satellite. It has an inclination of 
8.2598100e + 01°, a semi-major axis of 7.7788100e + 03 km, a B* value of 2.4820000e — 
04, and a period of 6.8277959^ + 03 seconds. It was propagated from 01-02-2014 to 01- 
09-2014 over 7 days. The least squares fit converged with an rms value of 2.11935 km. 
Therefore, it is concluded that the reason for the failure is the air drag effect. Figure 30 
represents the plots for position and velocity residuals over 7 days. 

SL-8 R/B, NORAD id 14966, has an inclination of 8.2972500e + 01°, a semi-major axis 
of 7.3620200e + 03 km, a B* value of 1.9141000e — 04, and a period of 6.2864617e + 03 
seconds. It was propagated from 01-01-2014 to 01-09-2014 over 9 days. The least squares 
fit converged with an rms value of 13.9352 km. Therefore, the air drag effect is the reason 
for the failure. 

Cosmos 2251, NORAD id 35981, has an inclination of 7.4059500e + 01°, a semi-major 
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Figure 30. The Residuals for Cosmos 2168, NORAD id 21782, over 9 Days 

axis of 6.9195900e + 03 km, a B* value of 2.9060000e — 03, and a period of 5.7283735e + 
03 seconds. It was propagated from 01-01-2013 to 01-14-2013 over 13 days. The least 
squares fit converged with an rms value of 86.645 km, which isn’t surprising because of its 
low altitude and high B* value. Therefore, the air drag effect is the reason for the failure. 
Figure 31 represents the plots for position and velocity residuals over 15 days. 

Cosmos 2297, NORAD id 23404, is in the close vicinity of the 14:1 resonance. It has 
an inclination of 7.1013200c + 01°, a semi-major axis of 7.2261100c + 03 km, a B* value 
of 2.4820000c — 04, and a period of 6.1131864c+ 03 seconds. It was propagated from 
01-01-2014 to 01-15-2014, 01-09-2014, and 01-07-2014 for three different time periods. 
All attempts failed because of the resonance. 

In conclusion, 1500 test cases, which fit the conditions mentioned in Section 3.1, were 
propagated using the low eccentricity KAM torus method. Then, each KAM torus propa¬ 
gation data was fitted to SGP4 and TLE prediction data associated with each test case over 
2 months. Initial success rate for convergence is 56.8667%. Then, 100 new test cases were 


Radial (U) 
Along-Track (V) 
Cross-Track (W) 


Radial (U) 

Along-Track (V) * 
Cross-Track (W) * 



94 






Position Residuals 


50 

0 

? 

~ ‘ 5 ° 

| -100 
IT 

-150 

-200 

23012 23014 23016 23018 23020 23022 23024 23026 

Time (Julian Dates) 


Radial (U) 
Along-Track (V) 
Cross-Track (W) 


0.05 

0 

-0.05 

-0.1 

-0.15 

-0.2 

23012 23014 23016 23018 23020 23022 23024 23026 

Time (Julian Dates) 


Velocity Residuals 


Radial (U) 
Along-Track (V) 
Cross-Track (W) 


Figure 31. The Residuals for Cosmos 2251, NORAD id 35981, over 15 Days 


selected from the failed test cases, and the same process was applied to 100 new test cases. 
The success rate of least squares fitting over 20 days instead of 2 months for the new test 
cases is 66%. Next, 34 failed test cases, which failed to converge for the second time, were 
individually analyzed. Seven test cases failed for the third time because their inclinations 
are in the close vicinity of polar inclination. An inclination of 86° caused the construction 
of the KAM torus to fail. The unstable forms of the Legendre polynomial recursions for 
the calculation of geopotential caused the high eccentricity problem. The solution is left 
for future studies. One test case failed to converge to the SGP4 and TLE data because of 
the inaccuracies in the TLEs. Ten test cases failed due to the air drag effect. Test cases with 
relatively high B* value and low altitude can only be propagated for short period of time 
because of effective air drag. Sixteen test cases failed to converge due to the resonance. 
There is no KAM torus near the resonance. Therefore, a clever approach should be taken 
to deal with the resonance, which is important for the generalization of the theory to all 
orbits. In addition, critical and polar inclination problems should also be addressed for the 
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generalization of the theory. 


4.2 Orbital Characteristics and Residuals 

The relationship between some orbital characteristics, such as semi-major axis, eccen¬ 
tricity, inclination, B* value, and period, and rms values yield valuable information about 
the theory. Because many factors degrade the accuracy and the dimensions that can be 
used to present the data are restricted to 3, straightforward presentation of data isn’t useful. 
Therefore, success ratio, which is successful attempts over the sum of the successful and 
failed attempts, for fixed time intervals versus orbital characteristics are used for the plots. 
The time intervals vary based on the orbital characteristics. However, there are 30 time 
intervals for all orbital characteristics. Each success ratio is an average value for a time in¬ 
terval, which is attributed to the starting point of the interval. This is important information 
for evaluating the plots including resonance data because the effects of the resonance don’t 
seem to be aligned. Therefore, additional explanations are provided for the plots. The data 
obtained from 1500 test cases were used for all plots in Section 4.1. 

The success ratio versus semi-major axis was plotted for LEO and GEO satellites be¬ 
cause there were few test cases at the MEO altitudes. For LEO satellites, the prime pertur¬ 
bation that degrades the accuracy for the theory is the air drag. The other parameter that 
degrades the accuracy for LEO satellites is the resonance. For GEO satellites, the prime 
perturbation is the third body mass. The other parameter that degrades the accuracy for 
GEO satellites is the 1:1 resonance. 

Figure 32 represents the success ratio versus semi-major axis for the test cases which 
have semi-major axes of 8000 km and smaller. The resonance points from left to right are 
the 29:2, 14:1, 27:2, 13:1, 37:3, and 12:1 resonances. Although 29:2 and 37:3 resonances 
do seem to be aligned with the decrease in the success ratio, the success ratio is an average 
over a time interval attributed to the starting point of the interval. For the 29:2 resonance 
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case, the success ratio average is from 7050 km to 7100 km, and the altitude, which is 
calculated from the center of Earth, for the 29:2 resonance is 7090.846 km for orbits with 
zero eccentricities. For the 37:3 resonance case, the success ratio average is from 7850 
km to 7900 km, and the altitude, which is calculated from the center of Earth, for the 37:3 
resonance is 7898.714 km for orbits with zero eccentricities. In addition, the 13:1 resonance 
doesn’t degrade the accuracy to zero for the interval between 7600 km and 7700 km. There 
are missing data for the interval between 7600 km and 7650 km. However, the resonance 
certainly degrades the accuracy. The success ratio increases with increasing altitude, which 
isn’t surprising, because air drag perturbation is more effective at low altitudes. 
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Figure 32. The Success Ratio versus Semi-Major Axis for LEO Objects over 2 Months 

The plot for the success ratio of GEO objects is surprising at the first glimpse because 
there are unexpected failures. Further examination yielded that three test cases created 
these failures. The author propagated these three test cases, which are Ius R/B, NORAD 
id 22316, Tvsat 2, NORAD id 20168, and familiar Galaxy 3, NORAD id 15308, from 
01-01-2014 to 03-01-2014. They converged with rms values of 13.3379 km, 10.1593 km, 
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and 14.4785 km, respectively, which proves that many failures can be related to the inac¬ 
curacies in TLEs, but this effort ignores this fact. Figure 33 represents the success ratio for 
test cases which have a semi-major axis between 42100 km and 42500 km before and after 
corrections are applied. The accuracy decreases rapidly after 42316 km. The convergence 
totally vanishes, or the accuracy becomes extremely low if the convergence is accomplished 
for 42416 km and bigger semi-major axis values due to the third body mass perturbations. 
The 1:1 resonance certainly degrades the accuracy. In addition, active payloads in geosyn¬ 
chronous orbit are maneuvered frequently to suppress the action of the resonance. 
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Figure 33. The Success Ratio versus Semi-Major Axis for GEO Objects 


Figure 34 represents the success ratio versus eccentricity for the test cases which have 
eccentricities between 0 and 0.01. In the plot, from 0 to 0.004 the accuracy is pretty lev¬ 
eled, except some ups and downs due to other orbital characteristics. However, after the 
threshold of 0.004 eccentricity is passed the success ratio decreases to approximately 30%. 
This proves how eccentricity is an important factor in terms of accuracy for the theory. 
Moreover, another analysis in Section 4.4 will examine the upper limit for the eccentricity 
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that can be modeled by the theory. 
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Figure 34. The Success Ratio versus Eccentricity over 2 Months 


Figure 35 represents the success ratio versus B * for the test cases which have B values 
between 0 and 0.01. In the plot, when the B* value of 0.003 is passed the success ratio 
for convergence decreases to approximately zero. The sharp decrease with increasing B* 
value shows the importance of the value for the theory. As mentioned previously, the B* 
value defines the susceptibility of near-Earth objects to air drag. The more effective the 
air drag is on a near-Earth object, the less accurate low eccentricity KAM torus theory 
becomes. The drastic effect of the B* value on the theory necessitates either an improved 
atmospheric model or the shorter propagation time interval. However, the first option may 
not introduce enough improvement in the accuracy due to the stochastic nature of air drag. 

One thousand five hundred test cases were pseudo-randomly selected under the conditions 
mentioned in details in Section 3.1, which include the close vicinity of polar and critical 
inclination, and eccentricity on the order of 10 3 and smaller. Figure 36 represents the 
success ratio versus inclination in order to check whether there are other inclination values 
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Figure 35. The Success Ratio versus Bstar over 2 Months 
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Figure 36. The Success Ratio versus Inclination over 2 Months 
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that lead to the resonance. In the plot, the inclinations of 20° and 80° seem to decrease the 
accuracy drastically, but this is an unexpected situation. Further examination yielded that 
three test cases for inclinations between 20° and 25°, two test cases for inclinations between 
25° and 30°, and three test cases for inclinations between 75° and 80° have created these 
failures. For inclinations 20°-25°, Delta 2 rocket body, NORAD id 21931, Pegasus rocket 
body, NORAD id 22491, and Scd 2, NORAD id 25504, have caused the failures. They all 
have low altitudes. Delta 2 R/B was propagated from 01-01-2013 to 01-10-2013, and the 
other two were propagated from 01-01-2014 to 01-10-2014. They converged to the SGP4 
and TLE data with rms values of 15.2938 km, 0.3447 km, and 0.2793 km, respectively. 
For inclinations 25°-30°, Tansei 1, NORAD id 4952, Delta 1 debris, NORAD id 10309, 
and Galex, NORAD id 27783, caused the failures. Delta 1 and Galex have low altitudes. 
Two of them were propagated from 01-01-2014 to 01-10-2014 with rms values of 269.943 
km, 0.2481 km, respectively. Although Delta 1 successfully converged, its semi-major 
axis, which is 6849.32, caused a big residual. Tansei 1 was also propagated from from 
01-01-2014 to 01-10-2014 in order to check whether attitude control, or station keeping 
maneuvers caused the failure. It was converged with an rms value of 0.5388 km. For 
inclinations 75°-80°, Scout X-4, NORAD id 871, and SL-21 rocket body, NORAD id 
29159, caused the failures. Two of them were propagated from 01-01-2014 to 01-10-2014 
with rms values of 0.9033 km, and 24.7307 km, respectively. Therefore, the polar and 
the critical inclinations are the only ones that fail the converge and degrade accuracy. The 
success ratio is quite leveled if the corrections are applied. Thus, there is no linear relation 
between the success ratio and the inclination. Moreover, another analysis in Section 4.4 
will examine the limits for the polar and critical inclinations that can be modeled by the 
theory. 

Figure 37 represents the success ratio versus the period for the test cases which have 
periods between 6000 seconds and 7160 seconds. The resonance points from left to right 
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are the 14:1, 13:1, and 12:1. The other order resonances, such as 27:2 and 37:3, can’t be 
detected in the plot. The two drastic decreases between 6400 seconds and 6600 seconds are 
due to the small number of test cases which all except one have high B* values. One of them 
failed to converge due to the inaccurate TLEs. For periods between 6440 seconds and 6480 
seconds. Delta 1 rocket body, NORAD id 21602, and Nimbus 2 debris, NORAD id 31579, 
reduced the accuracy. For periods between 6520 seconds and 6560 seconds, Cosmos 770, 
NORAD id 8325, which has a B* value of 2.845 x 10~ 3 , Meteor 3-2, NORAD is 19336, 
which has a B* value of 1.1514 x 10 -3 , Cosmos 1275 debris, NORAD id 14440, which has 
a B* value of 1.4665 x 10 2 , and SF-8 rocket body, NORAD id 11170, which failed due to 
the inaccurate TFEs because propagation from 01-01-2002 to 03-01-2002 yielded 16.7697 
km residual. Because there are few test cases between 6400 seconds and 6600 seconds, the 
resonance effect can’t be detected during this time period. Therefore, if the propagation 
time is decreased, the success ratio between 6400 seconds and 6600 seconds has a very 
similar trend as the success ratio between 6200 seconds and 6400 seconds, that is a linear 
increase with smaller variations. 

In conclusion, the resonance has emerged as a problem again in the broader picture. 
This section presents the analysis for 1500 test cases. For FEO objects, air drag degrades 
the success ratio of the least squares fit over 2 months to zero for altitudes below 422 
km. The accuracy also increases with increasing altitude. For GEO satellites, third body 
mass perturbation degrades the success ratio of the fit to zero for altitudes above 36038 
km. The only resonance that effects the accuracy is the 1:1 resonance for altitudes between 
35622 km and 36122 km. The eccentricity has appeared as an important parameter that 
degrades the accuracy. A sharp decrease in the accuracy occurs for eccentricities bigger 
than 0.004. The B* value has drastically decreased the accuracy of the theory. Moreover, 
all 16 test cases between the interval 0.0030 and 0.0033 have failed to converge. Therefore, 
the upper limit for B* for the theory can be defined as 0.003. However, the success rate for 
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Figure 37. The Success Ratio versus Period over 2 Months 

B* values between 0 and 0.001 can be easily enhanced by an improved atmospheric model. 
For the current theory, 217 out of 725 test cases have failed for the B* values between 0 
and 0.00033. The polar inclination and the critical inclination were proved to be the only 
problem that terminates the convergence of the least squares fit. Thus, there are no other 
inclination values that lead to failure of the convergence of the least squares fit. Moreover, 
there is no linear relation between the success ratio and the inclination. The success ratio 
based on the period is unsurprisingly very similar to the success ratio based on the altitude. 
There is a linear relation between the period and the success ratio of the convergence of the 
least squares fit. 

4.3 Some Samples from the Results 

The general performance of the theory has been presented so far. This section intro¬ 
duces some characteristic examples of the least squares fits for the LEO and GEO objects. 
Every test case has a plot that represents the accuracy of the low eccentricity KAM torus 
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theory assuming that SGP4 and TLE data are raw observational data. The reference frame 
for the plots is RAN frame of reference, see Section 2.4. However, this isn’t a fair evalu¬ 
ation of the accuracy of the theory because the basis frequencies and B * are continuously 
updated using SGP4 and TLE data. Section 4.5 introduces a fair comparison between 
SGP4 and the low eccentricity KAM torus method in terms of accuracy. Therefore, the 
plots introduced in this section present whether the numerically built torus can be fitted to 
observational data. If the torus is fitted to the observational data with high accuracy, the 
result is a compact numerical set of algorithms that are as accurate as numerical methods 
and as fast as general perturbation techniques, in terms of computational time. 
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Figure 38. The Residuals for Thorad Delta 1 Debris, NORAD id 8168, over 2 Months 


Figure 38 and Figure 39 represent the best 2 least squares fits for LEO objects. They 
are Thorad Delta 1 debris, NORAD is 8168, and Thorad Delta 1 debris, NORAD is 8140. 
The least squares fitting process yielded 3.6864262e — 01 km and 4.2330236e — 01 km, 
respectively. Although both of them have similar orbital characteristics, the difference 
in the accuracy is due to the B* values, which are 8.2174000<? — 05 and 9.3318000e — 
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Figure 39. The Residuals for Thorad Delta 1 Debris, NORAD id 8140, over 2 Months 


05, respectively. It should be also noted that Thorad Delta 1 debris, NORAD is 8140, is 
approximately at an altitude 25 km higher than Thorad Delta 1 debris, NORAD is 8168. 
Both plots present the bounded oscillatory behavior of the residuals over 2 months. They 
both have an altitude slightly above 1450 km. Thus, air drag isn’t that effective, and the 
quadratic structure of air drag usually appeared on along-track and radial components can’t 
be detected in the plots. 

Ops 6630 debris, NORAD id 9323, Cosmos 2251 debris, NORAD id 35469, and Tho¬ 
rad Agena D debris, NORAD id 4214, have least squares fits with residuals on the order 
of kilometers due to different reasons. Figure 40 represents the position and the velocity 
residuals of Ops 6630 debris for the least squares fit. Ops 6630 debris, NORAD id 9323, 
has an inclination of 9.7031200e + 01, a semi-major axis of 7.8056200e + 03, a B* value 
of 2.597 lOOOe — 04, and a period of 6.8631247e + 03 seconds. Its relatively high B* value 
leads to a residual of 5.02942 km over 2 months. However, the familiar quadratic structure 
of air drag that appears on the along-track and the radial components can’t be detected in 
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Figure 40. The Residuals for Ops 6630 Debris, NORAD id 9323, over 2 Months 


the plot because its B* value doesn’t exceed the limit that can be handled by the theory. 
However, Figure 41 represents the quadratic air drag effect on the least squares fit of Tho- 
rad Agena D debris. The along-track and the radial components seem to diverge from the 
desired oscillatory behavior of the residuals. Another factor that degrades the accuracy is 
the inaccuracies in the TLEs. Figure 42 represents the residuals of Cosmos 2251 debris for 
the least squares fit. In the plot, one of the TLEs is certainly inaccurate. 

GEO satellites don’t yield neat oscillatory behavior for the residuals because the third 
body mass perturbations aren’t added to the theory yet. However, third body mass pertur¬ 
bations can be calculated l ik e air drag perturbations, see Section 3.5. There will be at least 
one additional angle in the Fourier series to specify the motion of the Sun or the Moon. 
Three example test cases were chosen to present the characteristic fits for GEO satellites. 
Two of them have the best fits for GEO. The other one is one of the poorest fits. The poor 
fit proves the behavior of the residuals when the altitude of the satellite approaches to the 
threshold of 36000 km. Figure 43 represents the position and the velocity residuals of Ops 
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Figure 41. The Residuals for Thorad Agena D Debris, NORAD id 4214, over 2 Months 
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Figure 42. The Residuals for Cosmos 2251 Debris, NORAD id 35469, over 2 Months 
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Figure 43. The Residuals for Ops 9443 payload, NORAD id 11621, over 2 Months 


9443 payload, NORAD id 11621, for the least squares fit. The least squares fit for Ops 
9443 payload was converged with an rms value of 8.4179 km. Although this fit is one of 
the best fits for GEO satellites, its accuracy is far worse than the best fit for LEO objects 
because objects at high altitudes change their orbits towards the ecliptic plane of the solar 
system. Therefore, the third body mass perturbations should be added to the theory for bet¬ 
ter accuracy for GEO satellites. Figure 44 represents the position and the velocity residuals 
of Ops 9438 payload, NORAD id 10001, for the least squares fit. There are irregular pat¬ 
terns in the plot. However, third body mass perturbations are not the only source of error. 
Both station keeping maneuvers and inaccurate TLEs can cause these irregular residuals as 
well. Intelsat 4-F3 payload, NORAD id 5709, yielded one of the poorest fits. Intelsat 4-F3 
payload has a semi-major axis of 4.2349300e + 04 km and a period of 8.6732156e + 04 
seconds. It isn’t in the close vicinity of the 1:1 resonance. Thus, third body mass pertur¬ 
bations is the root cause of the poor fit. Figure 45 represents the position and the velocity 
residuals of Intelsat 4-F3 payload for the least squares fit. The propagation of Intelsat 4-F3 
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Figure 44. The Residuals for Ops 9438 payload, NORAD id 10001, over 2 Months 


4.4 Further Analysis for Orbital Characteristics 

The low eccentricity KAM theory is known to fail in the close vicinity of the polar and 
the critical inclinations and eccentricities above 10 3 . However, there is no experimental 
proof for the limits of the theory in terms of inclination and eccentricity. This section 
rigorously investigates those limits. For this section, 890 new test cases, which have an 
eccentricity of 10 3 and smaller, for polar inclination, 432 new test cases, which have an 
eccentricity of 10 3 and smaller, for the critical inclination, and 160 new test cases for 
eccentricity were retrieved from www. space-track. com servers, and propagated from 01- 
01-2014 to 03-01-2014. Another approach has been taken for the pseudo-random selection 
of the new test cases. For each case, the interval of interest was divided into 4 equal pieces. 
Equal number of random test cases were selected for each interval. 
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Figure 45. The Residuals for Intelsat 4-F3 payload, NORAD id 5709, over 2 Months 


Figure 46 represents the success ratio versus the eccentricity. The range for the eccen¬ 
tricity is between 0.01 and 0.05. This plot shows the upper limit of the eccentricity that can 
be modeled by the theory. The drastic decrease in the success ratio can be seen in the plot. 
The least squares fit clearly fails for the eccentricities of 0.0233 and bigger. The decrease 
in the success ratio for the eccentricity from 0 to 0.01 is 25%, see Section 4.2. However, 
the decrease in the success ratio for the eccentricity from 0.01 to 0.014 is approximately 
45%. 

Figure 47 represents the success ratio versus the inclination, which is in the close vicin¬ 
ity of the polar inclination. From 84° to 96° all 311 test cases have failed to converge. 
Therefore, the close vicinity of the polar inclination certainly terminates the fitting process 
completely. The current theory can’t model orbits which have inclinations between 84° 
and 96°. As mentioned previously, the solution to this problem is known, but it is left to 
future studies. This work is specifically beneficial for collision avoidance calculations and 
formation flight applications. 
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Figure 48 represents the success ratio versus the inclination, which is in the close vicin¬ 
ity of the critical inclination. The interval between 64.5° and 59.5° is where the success 
ratio becomes zero. Thus, the optimal region is the interval between 65.5° and 58°. Al¬ 
though it is hard to define the starting point of the optimal region due to the missing points, 
the upper limit is definitely 64.5°. 

The author also analyzed the 13:1 resonance. 156 tests cases failed between 6570 sec¬ 
onds and 6700 seconds. Delta 1 debris, NORAD id 6213, Jason debris, NORAD id 35414, 
and Cosmos 1691, NORAD id 35414, were the only survivors in the close vicinity of the 
13:1 resonance. Figure 49 represents the success ratio versus the period. This analysis 
provides valuable information for the behavior of the new theory near resonance. 

4.5 Comparison between KAM Torus Method and SGP4 

The least squares fitting process of the low eccentricity KAM torus to the SGP4 and 
TLE data requires continuous differential corrections of B* value, basis frequencies, and 
the the global variables. The author has written additional scripts to compare the new 
method to SGP4. First, a low KAM torus is built, and then, least squares is fitted to SGP4 
and TLE data over two months. Theory files, which are used to build the torus, and the 
basis frequencies are stored. Then, TLE of a future date, which is two months later than 
the starting date of the propagation, is propagated for an orbit by SGP4. This propagation 
is assumed to be the truth. Next, both the low eccentricity KAM torus and SGP4 are 
propagated over two months and residuals are calculated. The same process is repeated for 
the second time for other TLEs. 

Figure 50 represents the residuals over one period for the new theory. Figure 51 rep¬ 
resents the residuals over one period for SGP4. This comparison is made for the best 
least squares fitted near-Earth object, which has an rms value of 0.377006 km over two 
months. This test case is Thorad Delta 1 debris, NORAD id 8168. It has an inclination of 
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Figure 51. SGP4 prediction over 2 Months for the Best Fit Case 
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1.01334e-f-02°, a semi-major axis of 7.83195e + 03 km, and B* value of 8.21740^ — 05. 
The new theory yields approximately five times more accurate predictions than SGP4 does. 
The theory files for this comparison was created from 01-01-2013 to 03-01-2013. The date 
for the truth TLE is 05-01-2013. 


The KAM Theory Prediction For Ops 6630 



Figure 52. The Low Eccentricity KAM Torus Prediction over 2 Months for a Mean Case 

Figure 52 represents the residuals over one period for the new theory. Figure 53 repre¬ 
sents the residuals over one period for SGP4. This comparison is made for an average least 
squares fit, which has an rms value of 5.0332 km over two months. This test case is Ops 
6630 debris, NORAD id 9323. It has an inclination of 9.7031200e + 01°, a semi-major 
axis of 7.8056200e + 03 km, and B* value of 2.597lOOOe — 04. The new theory yields 
approximately three times more accurate predictions than SGP4 does. The theory files for 
this comparison was created from 01-01-2013 to 03-01-2013. The date for the truth TFE 
is 05-01-2013. 

The author also propagated the best fit case, NORAD id 8168, from 03-01-2013 to 05- 
01-2013 with approximately two days intervals. The result is very promising because there 
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Figure 53. SGP4 prediction over 2 Months for a Mean Case 
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is a linear error growth for one month with a small slope. Then, the error growth nearly 
doubles between 30 to 60 days. Figure 54 represents the residual growth versus period. 
The inaccuracies in SGP4 and TLE is also uncovered by Figure 54. 
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V. Conclusions and Recommendations 


This chapter presents the conclusions and the recommendations for future studies. The 
first section introduces the conclusions obtained from chapter 4. The limitations and appli¬ 
cability of the new theory is presented. For the second section, recommendations for future 
studies are given. 

5.1 Conclusions 

Although the title for this effort is “KAM Torus Orbit Prediction using SGP4 and TLE”, 
it has yielded more than only the possibility of orbit prediction combining the observational 
data and the theory. The new theory has proven to be a better substitute for SGP4. Its 
numerical set of algorithms provide higher accuracy and increase computational speed. 
The theory can be implemented in less than a minute on a home computer. The new theory 
has binary files associated with parameters, such as the periodic orbit, the geopotential and 
air drag perturbations, and second-order eccentricity perturbation. The new theory contains 
two parts. First, a data package which includes Fourier series. Second, a set of seven values, 
Q\ > Qh yi, y3, 34> ys, and B* at an epoch t 0 . Therefore, the new theory seems to be an 
alternate way to compress ephemerides. 

The low eccentricity KAM method appears to be promising for collision avoidance 
calculations because it yields accurate predictions. The linear residual growth of the new 
theory for short time intervals proves its applicability to collision avoidance calculations. 
Although it is hard to make a precise accuracy analysis without raw observational data, the 
error growth plot, see Section 4.5, proves that the low eccentricity KAM torus theory may 
rival numerical methods in terms of accuracy. In addition, the new theory can yield far 
better residuals with two simple modifications. First, the relatively inaccurate TLEs can be 
eliminated from the fitting process, see Figure 54. Second, the linear correction behavior 
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for the frequencies during the fitting process can be projected forward. Moreover, the new 
theory can reduce the dependence on NORAD. Because most ground stations depend on 
NORAD TLE, they must wait NORAD to publish the TLE. However, debris and most 
satellites can be propagated forward in time using TLE history with high accuracy using 
KAM torus theory. Formation flight is another possible application for the low eccentricity 
KAM torus theory. Since the accuracy of the theory increases for short time periods, the 
new method can be an onboard collision avoidance algorithm given that the satellite can 
receive position data. 

5.2 Future Studies 

There are some issues related to the new theory. As presented throughout chapter 4, 
resonance, high eccentricity, polar inclination, and critical inclination issues should be ad¬ 
dressed. This effort provides an extensive analysis of the new theory. The generalization 
of the theory to all orbits can be possible after these issues are solved. However, this effort 
can be used as future reference before applying theory on different debris and satellites. 
The solution to the the polar inclination problem is already known. The problem emerges 
due to the forms of the Legendre polynomial recursions which becomes numerically un¬ 
stable. The solution to this problem is provided in Wiesel’s text [64]. The atmospheric 
model should be improved. Russian GOST Atmosphere is a promising candidate because 
it requires less computation time and it is relatively accurate. It is built by empirical data 
obtained from Cosmos satellites. It includes solar flux and geomagnetic activity effects 
[60]. The application of the theory to orbits with high eccentricity and critical inclina¬ 
tion is an area of current interest. In addition, the actual accuracy analysis requires raw 
observational data, because SGP4 and TLE has intrinsic inaccuracies. 
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